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A 3-fluid hydrodynamic model for simulating relativistic heavy-ion collisions is introduced. Along¬ 
side with two baryon-rich fluids, the new model considers time-delayed evolution of a third, baryon- 
free (i.e. with zero net baryonic charge) fluid of newly produced particles. Its evolution is delayed 
due to a formation time r, during which the baryon-free fluid neither thermalizes nor interacts with 
the baryon-rich fluids. After the formation it starts to interact with the baryon-rich fluids and 
quickly gets thermalized. Within this model with pure hadronic equation of state, a systematic 
analysis of various observables at incident energies between few and about 160A GeV has been done 
as well as comparison with results of transport models. We have succeeded to reasonably reproduce 
a great body of experimental data in the incident energy range of I^iab — (1-160)A GeV. The list 
includes proton and pion rapidity distributions, proton transverse-mass spectra, rapidity distribu¬ 
tions of A and A hyperons, elliptic flow of protons and pious (with the exception of proton V 2 at 
40A GeV), multiplicities of pious, positive kaons, (j) mesons, hyperons and antihyperons, including 
multi-strange particles. This agreement is achieved on the expense of substantial enhancement of 
the interflow friction as compared to that estimated proceeding from hadronic free cross sections. 
However, we have also found out certain problems. The calculated yield of K~ is approximately 
by a factor of 1.5 higher than that in the experiment. We have also failed to describe directed 
transverse flow of protons and pion at i?iab > 40A GeV. This failure apparently indicates that the 
used EoS is too hard and thereby leaves room for a phase transition. 
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I. INTRODUCTION 

Relativistic nucleus-nucleus collisions at incident en¬ 
ergies Eiah ~ (10-40)A GeV presently attracts fecial 
attention, since the highest baryon densities [J, Q and 
highest relative strangeness 0 at moderate temperatures 
are expected in this incident energy range. The inter¬ 
est to this energy range has been recently also revived 
in connection with new experimental results from low- 
energy scanning SPS program 0 and the project of the 
new accelerator facility SIS300 at GSI 0- This domain 
of the nuclear phase diagram is less explored both ex¬ 
perimentally and theoretically as compared to the cor¬ 
responding extremes, i.e. cold compressed matter and 
baryon-free hot matter. Moreover, a critical point of the 
QCD phase diagram may occur to be accessible in these 
reaction s Ig , Q • There are already available experimen¬ 
tal data U pointing out that something interesting really 
happens in this energy range. In order to draw informa¬ 
tion on properties of hot and compressed nuclear matter 
from available experimental data, we need better under¬ 
standing of dynamics of nucleus-nucleus collisions under 
investigation. 
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A direct way to address thermodynamic properties of 
the matter produced in these reactions consists in appli¬ 
cation of hydrodynamic simulations to nuclear collisions. 
However, finite nuclear stopping power, revealing itself 
at high incident energies, makes the collision dynamics 
non-equilibrium and prevents us from application of con¬ 
ventional hydrodynamics especially at the initial stage of 
the reaction. Since the resulting non-equilibrium is quite 
strong, introduction of viscosity and thermal conductiv¬ 
ity does not help to overcome this difficulty, because by 
definition they are suitable for weak non-equilibrium. A 
possible way out is taking advantage of a multi-fluid ap¬ 
proximation to heavy-ion collisions, pioneered by Los- 
Alamos group [IJS and further developed at the Kur¬ 
chatov Institute ITlL , Frankfur t 11^ llH and 
GSI d d. The last development |lflL Il7|. i.e. the 
mean-field dynamics in the multi-fluid background, con¬ 
cerns mostly moderate energies around I GeV/nucleon 
and therefore will not be discussed here. 

The first extension of the original 2-fluid model H, , 
the two-fluid hydrodynamics with free-streaming radia¬ 
tion of pions, was advanced in ^3- The initial stage of 
heavy-ion collisions definitely is a highly non-equilibrium 
process. Within the hydrodynamic approach this non¬ 
equilibrium is simulated by means of a 2-fluid approx¬ 
imation, which takes care of the finite stopping power 
of nuclear matter [3113, and simultaneously describes 
the entropy generation at the initial stage. The radi¬ 
ated pions form a baryon-free matter in the mid-rapidity 
region, while two baryon-rich fluids simulate the propa- 
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gation of leading particles. The pions are the most abun¬ 
dant species of the baryon-free matter which may contain 
any hadronic, including baryon-antibaryon pairs, and/or 
quark-gluon species. 

First ap plic ations of the 2-fluid model with direct pion 
radiation (TH were quite successful for describing 
heavy-ion collisions in the wide range of incident ener¬ 
gies, from SIS to SPS. In these 3D hydrodynamic simu¬ 
lations the whole process of the reaction was described, 
i.e. the evolution from the formation of a hot and dense 
nuclear system to its subsequent decay. This is in dis¬ 
tinction to numerous other simulations, which treat only 
the expansion stage of a fireball formed in the course of 
the reaction, while the initial state of this dense and hot 
nuclear system is constructed from either kinetic simu¬ 
lations or more general albeit model-dependent as¬ 
sumptions (e.g. see 

However, the approximation of free-streaming pions, 
produced in the mid-rapidity region, was still irritat¬ 
ing from the theoretical point of view, in particular, be¬ 
cause the relative momenta of the produced pions and 
the leading baryons are in the range of the A resonance 
for the incident energies considered. This would imply 
that the interaction between the produced pions and the 
baryon-rich fluids should be strong. The free-streaming 
assumption relies on a long formation time of produced 
pions. Indeed, the proper time for the formation of the 
produced particles is commonly assumed to be of the 
order 1 fm/c in the comoving frame. Since the main 
part of the produced pions is quite relativistic at high 
incident energies, their formation time should be long 
enough in the reference frame of calculation to prevent 
them from interacting with the baryon-rich fluids. How¬ 
ever, this argument is qualitative rather than quantita¬ 
tive, and hence requires further verification. The first at¬ 
tempt to do this was undertaken by the Frankfurt group 
|13|. which started to explore an opposite extreme. They 
assumed that the produced pions immediately thermal- 
ize, forming a baryon-free fluid (or a “fireball” fluid, in 
terms of [H^I. and interact with the baryon-rich fluids. 
No formation time was allowed, and the strength of the 
corresponding interaction was guessed rather than mi¬ 
croscopically estimated. This opposite extreme, referred 
to as a (2-|-l)-fluid model and being not quite justified 
either, yielded results substantially different from those 
of the free-streaming approximation. This was one of 
the reasons why in subsequent applications the Frank¬ 
furt group neglected the interaction between baryon-free 
and baryon-rich fluids while keeping the produced pions 
thermalized [3113, thus effectively restoring the free- 
streaming approximation. However, the assumed imme¬ 
diate thermalization of the fireball fluid together with the 
lack of interaction with baryon-rich fluids still was not a 
consistent approximation. 

In the present paper we would like to introduce an ex¬ 
tension of the multi-fluid approach for simulating heavy- 
ion collisions, i.e. a 3-fluid model with formation time. 
This model is a straightforward extension of the 2-fluid 


model with radiation of direct pions [313 El and 
(2-|-l)-fiuid model [31313 We extend the above 
models in such a way that the created baryon-free fluid 
(which we call a “fireball” fluid, according to the Frank¬ 
furt group) is treated on equal footing with the baryon- 
rich ones. This implies that we allow a certain formation 
time for the fireball fluid, during which the matter of 
the fluid propagates without interactions. Moreover, we 
assume that the fireball matter gets quickly thermalized 
after its formation. The latter approximation is an en¬ 
forced one, since we deal with the hydrodynamics rather 
than with kinetics. The interaction between fireball and 
baryon-rich fluids is estimated based on elementary cross 
sections. A brief account of this model has been already 
reported in Ref. |^ . 


The formation time (r) is a conventional tool of the 
hadronic physics, which is associated with a finite time 
of string formation. It is incorporated in kinetic trans¬ 
port models such as UrQMD [^ and HSD |^. In dense 
medium the separate strings can already interact, this is 
the reason of introduction of junctions and formation of 
color ropes. In fact, this string interaction is a method 
of extending the treatment beyond the approximation of 
binary collisions which is inherent to the conventional 
Boltzmann equation. This interaction does not invali¬ 
date the concept of formation time but rather extends 
it to the case of multi-particle collisions. The concept of 
strings becomes irrelevant in the deconfined state of mat¬ 
ter. Here there are two points of view: (i) quarks and 
gluons are instantly produced accordingly to the QCD 
perturbation theor y, o r (ii) they are mediated by a coher¬ 
ent color field |3[3, i-®- the coherent color field 
is produced which subsequently decays into incoherent 
fluctuations - quarks and gluons. The first mechanism 
implies that r=0, while the second one still makes room 
for a finite r. Therefore, the fitted value of r may help to 
distinguish between these two mechanisms of the QGP 
production. 


The developed code allows calculations with various 
equations of states (EoS), which enter as a separate block 
of the code. We have started with the simplest, purely 
hadronic EoS which involves only a simple density de¬ 
pendent mean field providing saturation of cold nuclear 
matter at normal nuclear density no = 0.15 fm“^ and 
with the proper binding energy -16 MeV. This EoS is a 
natural reference point for any other more elaborate EoS. 
Much to our surprise, this trivial EoS turned out to be 
able to reasonably reproduce a great body of experimen¬ 
tal data. Taking advantage of the modern computers, 
substantial work has been also done on improvement of 
numerics of the model. These results are reported in this 
paper. 
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II. 3-FLUID HYDRODYNAMIC MODEL WITH 
DELAYED FORMATION 


The derivation of equations of the 3-fluid model pre¬ 
sented below is in fact only an illustration of physical 
assumptions which this model is based on. Indeed, the 
Boltzmann equation, from which this derivation starts, 
strictly speaking is inapplicable to a dense and strongly 
interacting matter. Therefore, the applicability of the 
3-fluid model to heavy-ion collisions is certainly an as¬ 
sumption which should be verihed in comparison with 
experimental data. However, this is a general situation 
with dense systems even in simpler cases. The basic rea¬ 
son for introduction of the 3-fluid approximaion is sim¬ 
ulation of the finite stopping power which is important 
at the formation stage of the initial hot and dense blob 
of nuclear matter. In this sense, it is an alternative to 
constructing this initial blob by means of either various 
kind kinetic transport models |2fl| or model assumptions 

mill. 


A. Basic Formulation 

Unlike the conventional hydrodynamics, where local 
instantaneous stopping of projectile and target matter 
is assumed, a specific feature of the dynamic 3-fluid de¬ 
scription is a finite stopping power resulting in a counter¬ 
streaming regime of leading baryon-rich matter. Exper¬ 
imental rapidity distributions in nucleus-nucleus colli¬ 
sions support this counter-streaming behavior, which can 
be observed for incident energies between few and 200A 
GeV. The basic idea of a 3-fluid approximation to heavy- 
ion collisions caiii is that at each space-time point 
X = (t, x) the generally nonequilibrium distribution func¬ 
tion of baryon-rich matter, /br(a;,p), can be represented 
as a sum of two distinct contributions 


where Ca denote collision terms between the constituents 
of the three fluids. We have omitted intra-fluid collision 
terms, like Cp{fp, /p), since below they will be canceled 
anyway. The displayed inter-fluid collision terms have 
a clear physical meaning: Cp/t{fp, ft), Cp/t{fp/t, f{), and 
Cf(/p/t, /f) give rise to friction between p-, t- and f-fluids, 
and the term C'f(/p, ft) takes care of particle production 
in the mid-rapidity region. Note that up to now we have 
done no approximation, except for hiding intra-fluid col¬ 
lision terms. 

Let us proceed to approximations which justify the 
term “fluids” having been used already. We assume that 
constituents within each fluid are locally equilibrated, 
both thermodynamically and chemically, i.e. that fa 
are equilibrium distributions. In particular, this implies 
that the intra-fluid collision terms are indeed zero. This 
assumption relies on the fact that intra-fluid collisions 
are much more efficient in driving a system to equilib¬ 
rium than inter-fluid interactions. As applied to the hre- 
ball fluid, this assumption requires some additional com¬ 
ments, related to the concept of a finite formation time. 
During the proper formation time r after production, the 
fireball fluid propagates freely, interacting neither with it¬ 
self nor with the baryon-rich fluids. After this time inter¬ 
val, the fireball matter starts to interact with both itself 
and the baryon-rich fluids and, as a result, locally ther- 
malizes. Being heated up, these three fluids may contain 
not only hadronic and but also deconfined quark-gluon 
species, depending on the EoS used. 

The above assumption suggests that interaction be¬ 
tween different fluids should be treated dynamically. To 
obtain the required dynamic equations, we first integrate 
the kinetic Eqs. over momentum and sum over 

particle species with weight of baryon charge. This way 
we arrive to equations of the baryon charge conservation 

dMx) = 0, (5) 


fhv{x,p) = fp{x,p) + ft{x,p), (I) 

initially associated with constituent nucleons of the pro¬ 
jectile (p) and target (t) nuclei. In addition, newly pro¬ 
duced particles, populating the mid-rapidity region, are 
associated with a fireball (f) fluid described by the distri¬ 
bution function f[(x,p). Therefore, the 3- fluid approxi¬ 
mation is a minimal way to simulate the finite stopping 
power at high incident energies. Note that both the 
baryon-rich and fireball fluids may consist of any type 
of hadrons and/or partons (quarks and gluons), rather 
than only nucleons and pions. However, here and below 
we suppress the species label at the distribution functions 
for the sake of transparency of the equations. 

With the above-introduced distribution functions fa 
(a =p, t, f), the coupled set of relativistic Boltzmann 
equations looks as follows: 

= Cp(/p,/t) + Cp(/p,/f), (2) 

Pf,d>ift{x,p) = Ct(/p,/t) + a(/t,/f), (3) 

P^d>fMx,p) = Cf(/p, ft) + Ciifp, ft) + Ctift, /f),(4) 


for a =p and t, where = UaU^ is the baryon current 
defined in terms of baryon density Ua and hydrodynamic 
4-velocity normalized as Ua^iU^ = 1. Eq. @ implies 
that there is no baryon-charge exchange between p- and 
t-fluids, as well as that the baryon current of the hre- 
ball fluid is identically zero, = 0. Integrating kinetic 
Eqs. @-0 over momentum with weight of 4-momentum 
p'^ and summing over all particle species, we arrive at 
equations of the energy-momentum exchange for energy- 
momentum tensors Tff' of the fluids 

d^T^''{x) = -F;{x)+Ff^{x), (6) 

d^^Trix) = -Ff{x)+Ff,{x), (7) 

d^T^\x) = F;{x) + Ff(x) - Ff^{x) - Fff{x), (8) 

where the F^ are friction forces originating from inter¬ 
fluid collision terms in the kinetic Eqs. @-0. Fp and 
Ff in Eqs. 0-0 describe energy-momentum loss of 
baryon-rich fluids due to their mutual friction. A part of 
this loss \Fp — Ff\ is transformed into thermal excitation 
of these fluids, while another part {F^ + Ff) gives rise 


4 


to particle production into the fireball fluid (see Eq. (O). 
Ffp and are associated with friction of the fireball fluid 
with the p- and t-fluids, respectively. Note that Eqs. m- 
© satisfy the total energy-momentum conservation 

= 0. (9) 

As described above, the energy-momentum tensors of 
the baryon-rich fluids {a =p and t) take the conventional 
hydrodynamic form 


Tr = (£„ + p„)< < - 


( 10 ) 


in terms of the proper energy density, and pres¬ 
sure, Pq,. For the fireball, however, only the energy- 
momentum tensor of the formed matter is of practical in¬ 
terest for us. Since we treat the fireball matter as a fluid, 
we have nothing to do but assume that it is formed al¬ 
ready thermalized. This assumption implies that its ther- 
malization time is essentially shorter than its formation 
time. Moreover, this assumption is in the spirit of other 
assumptions made: we distinguish the main nonequilib¬ 
rium associated with finite stopping power and consider 
all the rest of dynamics within local equilibrium. Thus, 
only the formed (and by assumption thermalized) part 
of the energy-momentum tensor is described by this hy¬ 
drodynamic form 

+ Pf) (11) 


Its evolution is defined by an Euler equation with a re¬ 
tarded source term 

= -Pf'^(x) - Pf'((x) 

-I- J d'^x'S'^ (x - x' - Uf{x')t) [Fp{x') + {x')] ,{12) 

where r is the formation time, and 


U^ix') 


Up{x')+uUx') 

\Up{x') + Ut{x')\ 


(13) 


The above discussed free-propagating fireball mat¬ 
ter is, of course, a certain approximation for treating 
free-streaming particles, which this matter consists of. 
These produced free-streaming particles are character¬ 
ized by a distribution function, /(p), in the momen¬ 
tum space, which is determined by their production 
cross sections, e.g. p^dufsiN^-nX/d?p. In particular, this 
means that particles get formed at different time in¬ 
stants TparticieTparticie, depending on the particle velocity 
(7particie IS the gamma factor of the particle). This par¬ 
ticle formation is governed by a particle formation time 
Tparticie which differs from our r we use in the hydrody¬ 
namic formulation. Our r has a meaning of the parti¬ 
cle formation time Tparticie averaged over the distribution 

f{p)- 


r j r d^p 

X "f{ = J "Tparticie 7particle/(p) / J / (p) ) (14) 


where 7 f is the gamma factor of the unformed fireball 
matter. Eq. da shows that r is always longer than 
Tparticie- In particular, if the fireball matter consists of pi- 
ons, which have approximately thermal distribution with 
temperature « 100 MeV, then r « 2Tparticie. 

The free-streaming particles are not only formed at 
different time instants but also in different space points 
because of the same /(p) distribution. Thus, another 
approximation we have done consists in neglecting this 
spacial spread. We could take it into account by chang¬ 
ing the ^-function in the r.h.s. of Eq. 03 to a smooth 
distribution of formed particles. This would result in spa¬ 
tial smearing of formation of the fireball fluid. However, 
at the same time it would make the model formulation 
more complex, since this smooth distribution of formed 
particles should depend on the initial /(p) distribution 
and, hence, on the cross sections, like p'^duNN^ T,x/d?p. 
Moreover, it would be another source of uncertainty, since 
the NN fX process is not the only relevant one. 
Therefore, for the first estimate of effects of the formation 
time we avoid these complications. 

The residual part of (the free-propagating one) is 
defined as 


is a free-propagating 4-velocity of the produced fireball 
matter, which is evidently a time-like 4-vector. In fact, 
this is the velocity at the moment of production of the 
fireball matter. According to Eq. 03, the energy and 
momentum of this matter appear as a source in the Euler 
equation only later, i.e. after the time interval U^t upon 
the production, and in different space point x'—Uir'(a;') r, 
as compared to the production point x'. From the first 
glance, one can immediately simplify the r.h.s. of Eq. 
d by performing integration with the i5-function. How¬ 
ever, this integration is not so straightforward, since the 
expression under the i5-function, x — x' — Uf{x')t = 0, 
may have more than one solution with respect to x' . The 
latter would mean that the matter produced in several 
different space-time points x' is simultaneously formed in 
the same space-time point x. This is possible due to the 
nonlinearity of the hydrodynamic equations. 


The equation for can be easily obtained by tak¬ 

ing the difference between Eqs. © and (I12II . If all the 
fireball matter turns out to be formed before freeze-out, 
then this equation is not needed. Thus, the 3-fluid model 
introduced here contains both the original 2-fluid model 
with pion radiation [3113 and the (2-|-l)-fluid model 
[l^ 13 . I3| as limiting cases for r ^ oo and r = 0, 
respectively. 

B. Friction between Baryon-Rich Fluids and their 
Unification 

The nucleon-nucleon cross sections at high energies 
are strongly forward-backward peaked. This fact, which 
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originally served as justification for subdividing baryonic 
matter into target and projectile fluids, was used in la 
to estimate the friction forces, and proceeding 
from only NN elastic scattering. Later these friction 
forces were calculated based on (both elastic and in¬ 
elastic) experimental inclusive proton-proton cross sec¬ 
tions. In the present calculations we use the following 
form of the projectile-target friction 

K = [« - Dp + (wK + <) De] , (16) 

a =p or t, p =t and t =p. Here, p% denotes a kind 
of ’’scalar densities” of the p- and t-fluids (see details 
below), 

Dp/E = '^N (Tp/£;(Spt), (17) 

where rriAr is the nucleon mass, Spt = (up -I- is 
the mean invariant energy squared of two colliding nu¬ 
cleons from the p- and t-fluids, 

Kei = [spt(spt - 4m^)]i/V2m^ (18) 


because it does not take into account : 

(i) various mesonic (and maybe, quark and gluonic) 
species produced in the collision 

(ii) possible multiparticle interactions which are quite 
probable in the dense medium, 

(iii) possible medium modifications of cross sections and 
effective masses, and 

(iv) quark and gluon interactions, if the system happens 
to undergo the phase transition to the quark-gluon phase. 

(v) Moreover, even experimental cross sections between 
hadronic species are only poorly known. As for the non- 
perturbative quark-gluon phase, this is the range of pure 
speculations. 

In view of theses uncertainties, it is reasonable to make 
provision for tuning the above friction. For this purpose 
we introduced tuning factors ^(spt) in the scalar densities 
of the p- and t-fluids 

pii^pt) = (pr-+ a(spt) 

+ 2 (Pa + Pa) Cg('Spt)j (21) 


is the mean relative velocity of the p- and t-fluids, and 
o’p/£;(spt) are determined in terms of nucleon-nucleon 
cross sections integrated with certain weights (see E 
mil for details): 


o'p(spt)= / daNN^NX — cos 9 c 

•l9cm<7r/2 V 

0'EiSpt)= f daNN^Nxi^ -) • 

■lecm<7r/2 V An J 


Pont N 

T I 


(19) 

( 20 ) 


Here the integration is restricted to the forward hemi¬ 
sphere (0cm < 7’‘/2) of the center-of-mass scattering an¬ 
gles 0cm, Pin = (spt/4-m^)i/2 and An = Sp{^/2 are the 
in-coming momentum and energy of the nucleon before 
the scattering in the NN c.m. frame, respectively, and 
Pout and Aut are the corresponding out-coming quan¬ 
tities. (Tp(spt) is a kind of a transport cross section, 
which is nonzero at any physical Spt, as it is seen from 
Eq. (I19II . At the same time, the (T£;(spt) quantity, which 
is responsible for the fireball production, is zero for Spt 
below the inelastic threshold. In particular, the latter 
feature makes Up = 0, and hence the friction A = 0, 
at vanishing relative velocity Up — Ut = 0. The overall 
factor in Eq. Util) is associated with unification of p- 
and t-fluid into a single one (see Eq. (I^ l. when their 
relative velocity gets small enough. 

The above friction fTC|) is a certain extension of that 
derived in and used in m HI HI HI in This 
original derivation was performed under assumption that 
baryon-rich fluids consist of only nucleons, or baryons, as¬ 
suming that baryon-baryon cross sections are similar to 
proton- pro ton ones. Therefore, in the original expression 
for A 0, 0, 0,1^ modified scalar densities 

of the p- and t-fluids are substituted by corresponding 
baryon densities (with A = Cg = 1 and = 1, see be¬ 
low). However, this original calculation is incomplete, 


where p^^', Pa and p® are scalar densities of all 

baryons, all mesons, quarks and gluons, respectively, de¬ 
fined in the conventional way, i.e. 


Pa{x) = ma 


d^p 

Po 


/a(a:,p) 


( 22 ) 


for the a species of mass nia which are pa{T(x), p,b{x)) in 
the thermodynamic limit. These scalar densities (but of 
course, only nucleon ones) result from the original deriva¬ 
tion of Ref. [^. We just extended the recipe of Ref. 
to other species. In fact, in the nonrelativistic case the 
scalar density is identical to the usual particle-number 
density, since the mass ma in Eq. (I22II is canceled by 
Po in the denominator. Moreover, in the relativistic case 
the scalar density is a natural result, since it is a Lorentz 
invariant contrary to the particle-number density, which 
is the 4th component of a 4-vector. However, in the rela¬ 
tivistic case we run into trouble with gluons and quarks. 
If one uses current masses for gluons and quarks, the cor¬ 
responding scalar densities are either identically zero (for 
gluons) or negligibly small (for u and d quarks). To over¬ 
come this problem we imply that either thermal effective 
masses or constituent masses, depending on particular 
model, should be used in the scalar densities, or they 
should be calculated beyond the quasiparticle approxi¬ 
mation (e.g. on the lattice). 

Factors like 2/3 and 1/3 in Eq. 1211) take into account 
the assumed scaling of cross sections in accordance with 
the naive valence-quark counting. Note that the original 
free cross section in oa was estimated for the proton- 
proton pair. In Eq. m different tuning factors are 
introduced for hadronic and quark-gluon phase: A and 
respectively. If the system occurs in pure hadronic 
phase, directly scales the cross sections op/E- If the 
system happens to be in a mixed phase, like in the cross¬ 
over phase transition, both A and participate in the 
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cross-section scaling. Note that originally the scalar den¬ 
sities depend only on x, cf. Eq. The introduced 

s-dependent tuning factors make the tuned scalar den¬ 
sities also s-dependent. In general, and could 
be also assumed to be temperature dependent, since the 
ap/E are obtained by averaging over certain thermal dis¬ 
tributions. Nevertheless, we avoid to do this in order to 
keep the number of “fitting degrees of freedom” as low 
as possible. 

When the baryon-rich fluids get decelerated enough, it 
is reasonable to unify them into a single fluid. This limit 
is not automatically present in the m odel formulation and 
therefore should be stiDulated |lO(i| . For this purpose we 
introduce an auxiliary function 

^ = l_exp[-(V;Pf/Al/)4] (23) 

where is the mean relative velocity of the p- and t- 
fluids of Eq. (I18II , and AV is the characteristic velocity 
of particles inside the flow (the Fermi velocity of nucleons 
in cold matter or the thermal velocity of particles in hot 
matter) 




1 = max 

p,t 


ep 3r 1 
mjv ’ 2mN j ’ 


(24) 


where ep is the Fermi energy corresponding to the baryon 
density of the fluid, T is temperature of the fluid, and the 
maximum value is searched over quantities correspond¬ 
ing to the two overlapped (p- and t-) fluids. Thus, at 
AF, -i? = 1 and the purely two-fluid regime is 
realized, while at <C AF, 19 = 0 and the one-fluid 
regime takes place, which implies complete equilibration 
of the overlapping fluids. During the numeric simulation 
the two-fluid and one-fluid solutions are mixed in propor¬ 
tion i9 and 1 — "d, respectively, thus providing a smooth 
approach to the one-fluid limit. This interpolation pro¬ 
cedure concerns the pressure Pa and hydrodynamic 3- 
velocity Vq,, i.e. those quantities that are required for 
the hydrodynamic transport and are non-additive: 

ji 

Pa=l9Pa + (l-l9)—Ptot, (25) 

‘^tot 

V„ = -h (1 -'(?)vtof (26) 


Here Pa and Vq, are interpolated quantities for the a 
fluid, and Ptot, Vtot, n-tot, etc. are quantities derived 
from the total baryon density = Jp -I- and the Op 
components of the total baryon-rich energy-momentum 
tensor -|- assuming that p- and t-fluids 

are unified. The difference between Pa and Ptot, and be¬ 
tween Vq, and Vtot, respectively, is as follows. From the 
solution of hydrodynamic equations ®-0 and m at 
each time step we obtain hydrodynamic quantities Jq, 
and same time, we need, in partic¬ 

ular, the pressure and hydrodynamic 4-velocity in order 
to proceed to the next time step. We have to calculate 
them based on the above hydrodynamic quantities and 


the EoS. Here we can proceed in two ways. The con¬ 
ventional way of the multi-fluid hydrodynamics consists 
in calculating the pressure and 4-velocity for each fluid, 
i.e. solely based on the Jq and quantities related to 
this fluid. This way we arrive at Pq and Vq, which com¬ 
pletely preserve the multi-fluid character of the solution. 
Alternatively, we can completely abandon the multi-fluid 
nature and assume that we deal with a locally equilib¬ 
rium piece of baryon-rich matter, which is characterized 
by J^at and and then calculate the pressure and 4- 
velocity of this “unified” piece. Thereby we arrive at Ptot 
and Vtot, which describe unified baryon-rich fluids as if 
they are mutually stopped. 

This unification procedure was first proposed by the 
Los-Alamos group Q and then used in subsequent ap¬ 
plications of the multi-fluid dynamics Pin in m. 
Evidently, this procedure also affects the stopping power 
and therefore can be added to the above list of uncertain¬ 
ties associated with the friction force Pq. Other criteria 
of unification may produce different observable stopping, 
as it was shown in Ref. |15|. However, we prefer to keep 
the original criterion, keeping in mind that this is an im¬ 
portant part of the stopping prescription. 

As compared to the already conventional unification 
procedure described above, we have also introduced i9^ 
factor in the friction force itself, cf. Eq. It is justi¬ 

fied, since only the i9 fraction of each baryon-rich fluid re¬ 
mains in the two-fluid regime, as it follows from the above 
discussion. Therefore, in order to be consistent with the 
above unification procedure, we should keep only this i9 
fraction of the density of each fluid, in the friction 
term. 

As it has been already mentioned, the way, in which 
the unification procedure is realized, affects the stopping 
of the nuclear matter. In fact, it is possible to avoid this 
artificial unification procedure. As it was shown in Refs. 
Elli, the formulation of the multi-fluid dynamics in 
terms of mean fields rather than of the EoS results in 
automatic unification of mutually stopped matter. How¬ 
ever, the problem of such mean-field formulation is that 
the mean fields become unrealistically strong at high rela¬ 
tive velocities of counter-streaming matter, whereas they 
should die out due to momentum dependence of self en¬ 
ergies, as it was advocated in the HSD model Ref. [^ . 
Therefore, for the domain of high incident energies the 
present approach, based on a separate EoS in each fluid 
and complemented by the unification procedure, is cer¬ 
tainly preferable as compared to the pure mean-field for¬ 
mulation of Refs. |ial3- While the latter formulation 
is definitely advantageous at moderate incident energies 
of the order of lA GeV. 

Eqs. 0-0 and lUTl) . supplemented by a certain EoS 
and expressions for friction forces F'^, form a full set of 
equations of the relativistic 3-fluid hydrodynamic model. 
The only quantity, we still need to define in terms of 
hydrodynamic variables and some cross sections, is the 
friction of the fireball fluid with the p- and t-fluids, 
and 
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C. Interaction between Fireball and Baryon-Rich 
Fluids 

Our aim here is to estimate the scale of the friction 
force between the fireball and baryon-rich fluids, similar 
to that done before for baryon-rich fluids To this 

end, we consider a simplified system, where all baryon- 
rich fluids consist only of nucleons, as the most abundant 
component of these fluids, and the fireball fluid contains 
only pions. 

For incident energies from lOA (AGS) to 200A GeV 
(SPS), the relative nucleon-pion energies are in the res¬ 
onance range dominated by the A-resonance. To esti¬ 
mate this relative energy we consider a produced pion, 
being at rest in the center of mass of the colliding nu¬ 
clei, p = 0, 0,0}cm- Baryon-rich fluids decelerate 

each other during their inter-penetration. This means 
that the nucleon momentum q should be smaller than 
the incident momentum, go = < rnN^cm, where 7cm 

is the gamma factor of the incident nucleon in the c.m. 
frame. Calculating the invariant relative energy squared 
s = (p + at Aiab = 158A GeV, we obtain <1.8 
GeV. This range of s precisely covers the resonance re¬ 
gion, 1.1 GeV < < X g Qgv [13. At Fiiab = lOA 

GeV we arrive at < 1.3 GeV, which is also within 

the resonance region. At even lower incident energies the 
strength of the fireball fluid becomes so insignificant, as 
compared with thermal mesons in the p- and t-fluids, 
that the way of treatment of its interaction with the 
baryon-rich fluids does not essentially affect the observ¬ 
ables. For the same reason we do not apply any special 
prescription for the unification of the fireball fluid with 
the baryon-rich fluids, since this may happen only at rel¬ 
atively low incident energies Aiab < lOA GeV. 

The resonance-dominated interaction implies that the 
essential process is absorption of a fireball pion by a p- or 
t-fluid nucleon with formation of an i?-resonance (most 
probably A). This produced i?-resonance still belongs to 
the original p- or t-fluid, since its recoil due to absorption 
of a light pion is small. Subsequently this i?-resonance 
decays into a nucleon and a pion already belonging to 
the original p- or t-fluid. Symbolically, this mechanism 
can be expressed as 

N°‘ + ^ R°‘ ^ N°‘ + ■K°‘. 

As a consequence, only the loss term contributes to the 
kinetic equation for the fireball fluid. 

Proceeding from the above consideration, we write 
down the collision term between fireball-fluid pions and 
a-fluid nucleons {a =p or t) as follows 

QiUJf) = - / ^ ft'^\p)Uq), (27) 

J qo 


is the rate to produce a baryon i?-resonance, and 
(s) is the parameterization of experimental pion- 
nucleon cross-sections 113. Here, only the distribution 
function of formed (and hence thermalized) fireball pi¬ 
ons, enters the collision term, since the non-formed 

particles do not participate in the interaction by assump¬ 
tion. 

Multiplying C{{fa^ fi) by the 4-momentum p'' and in¬ 
tegrating the result over momentum, we arrive at 


FtUx) = f Uq) 

qo Po 


W 


Nn^R( 


(sfa) f f (fq 


niTrUf 


qo 


U{q) 


(fp 

Po 


0 1/ f(eq) 

p p ft 


= Dfa 




(eq)Oi/ 


ip) 

(28) 


where we substituted p^ and s by their mean values, 
(pP) = and Sfa = {m.^Ui + mNUa)'^, and introduced 
the transport coefficient 


Dta = 


W 


Nn^R 


(sta) 


= K 


rriN'm-K 


. fa. Ntt^R 
rel ^tot 


(sfc). (29) 


Here, 17^7 = [(s^ - rn% - m^)^ _ {2mNm^) 

denotes the mean invariant relative velocity between 
the fireball and the a-fluids. Thus, we have expressed 
the friction in terms of the fireball-fluid energy- 
momentum density (of only pions as yet), the scalar 
density pa of the a-fluid (of only nucleons as yet), and a 
transport coefficient Dfa- Note that this friction is zero 
until the fireball pions are formed, since = 0 dur¬ 

ing the formation time r. 

In fact, the above treatment is an estimate of the fric¬ 
tion terms rather than their strict derivation. All the un¬ 
certainties mentioned in the previous subsection are well 
applied to the case under consideration. This peculiar 
way of evaluation is motivated by the form of the final 
result (1281 . An advantage of this form is that m,r and any 
other mass do not appear explicitly, and hence it allows 
a natural extension to any content of the fluid, including 
deconfined quarks and gluons, assuming that Dta repre¬ 
sents just a scale of the transport coefficient. Performing 
such extension, we assume that represents the total 
energy-momentum density of the fireball-fluid and pa is 
the total scalar density of the a-fluid, i.e. pa = in 
terms of Eq. m- Here we have omitted tuning ^ factors 
in the scalar density pa, since in this case their effect is 
very similar to that of the formation time, which switches 
on/off the interaction at various stages and thereby ef¬ 
fectively changes its strength. 


D. Freeze-Out 


where s = {p + g)^, 

= ^\/is-m%- mlY - (s) 


The hydrodynamic simulation is terminated by a 
freeze-out procedure. Though this method (as applied to 
high-energy physics) was first proposed almost 50 years 






ago [^, this is still an actual problem which is actively 
discussed. The method was intuitively clear and easily 
applicable. However, as was shown by Cooper and Frye 
|33 , the Milekhin’s method violates the energy conserva¬ 
tion. To remedy the situation, they proposed their own 
recipe. The Cooper-Frye recipe was not free of prob¬ 
lems as well. It gives negative contribution to the particle 
spectrum in some kinematic regions in which the normal 
vector to the freeze-out hyper-surface is space-like. This 
negative contribution corresponds to frozen out particles 
returning to the hydro phase. Cut off of this negative con¬ 
tribution again returns us to the violation of the energy 
conservation. The get rid of this negative spectrum, there 
was proposed a modification of the Cooper-Frye recipe 
based on a cut-Jiittner distribution [^13 EH 113 ■ In 
this distribution the part of the Jiittner distribution that 
gave the negative spectrum is simply cut off. To preserve 
the particle and energy conservation, the rest of Jiittner 
distribution is renormalized, effectively resulting in a new 
temperature and chemical potential (so called ” freeze-out 
shock”). In fact, this cut-Jiittner recipe has no physi¬ 
cal justification, except for practical utility. Moreover, 
the cut-Jiittner recipe is not supported by schematic ki¬ 
netic treatment of the transitional region from hydro 
regime to that of dilute gas and looks like a violence to the 
nature, making the freeze-out procedure uncontrollable. 
Recently there was proposed a new freeze-out recipe, a 
cancelling-Jiittner distribution [^, which complies with 
results of schematic kinetic treatment [^. However, very 
recently the authors reported that this cancelling- 
Jiittner distribution is satisfactory only for the space¬ 
like freeze-out, while it fails for the time-like one. It 
should be stressed that this was precisely the schematic 
kinetic treatment. This region, where the transition from 
highly collisional dynamics to the collisionless one occurs, 
is highly difficult for the kinetic treatment and hardly al¬ 
lows any justified simplifications. 

All the above considerations of the freeze-out process 
proceeded from assumption of existing some continu¬ 
ous hyper-surface separating the hydro system from the 
frozen-out gas. Conservation conditions on such hyper¬ 
surface are constructed in analogy with shock front in 
hydrodynamics. From the practical point of view, it 
would mean that we should first run the hydro calcu¬ 
lation without any freeze-out and only after that look for 
a hyper-surface, where the freeze-out criterion is met. In 
practice our hydro simulation proceeds in different way. 
The freeze-out criterion is checked continuously during 
the simulation. If some parts of the hydro system meet 
this criterion, they decouple from the hydro calculation. 
The frozen-out matter escapes from the system, removing 
all the energy and momentum accumulated in this mat¬ 
ter. Therefore, it pro duce s no recoil to the rest of still 
hydrodynamic svstem [l07| . In particular, it means that 
the boundary condition on the free surface (between the 
hydro system and vacuum) is not kept at the same po¬ 
sition as in the calculation without freeze-out but moves 
inside the system. It affects not only the system surface 


but also the interior. The freeze-out process looks like an 
evaporation (or fragmentation, on account of final-size 
grid) first from the system surface and then as a volume 
fragmentation of the system residue. 

In view of above said, we prefer to avoid not quite 
justified complications of the freeze-out procedure and 
make use of the simplest (how ever equally unjustified) 
original choice of Milekhin with corrected treatment 
of the energy conservation. The freeze-out criterion we 
use is 

^tot ^ ^frz; (^ 0 ) 

where 

£tot = (31) 

\ / proper 

is the energy density of all three fluids in the proper ref¬ 
erence frame, where all nondiagonal components of the 
total energy-momentum tensor are zero, 

(t^'' + rr + = 0, (32) 

V / proper 

and Efrz is the critical freeze-out energy density. If the 
freeze-out criterion is met in some space-time point, all 
three fluids in this point get frozen out and removed from 
the hydrodynamic evolution. In fact, we freeze out the 
fluids in tiny portions, i.e. droplets. This is allowed and 
even implied by the numerical scheme we use (see App. 
El). Each droplet gets frozen out in its proper reference 
frame. In terms of the hypersurface, our freeze-out hyper¬ 
surface is discontinuous and consists of small fragments 
with normal vectors of local 4-velocities. Any transport 
is prohibited on the non-existent parts (parallel to the 
flow velocity) of the surface. In this sense, our freeze- 
out is more similar to a continuous fragmentation of the 
system rather than to occurrence of the shock front. 

A frozen-out droplet of the a-fluid (let it be marked 
as ia) is still characterized by some temperature, baryon 
and strange chemical potentials corresponding to a non¬ 
gas EoS (involving some mean fields) used in the hydro 
calculation. This is not suitable for calculation of the 
spectrum of observable particles. First we should re¬ 
lease the energy stored in mean fields. To do this, we 
recalculate temperature (T®“(gas)j^ baryon and 

strange chemical potentials corresponding to 

the hadronic gas EoS proceeding from conservations of 
total energy-momentum, baryon and strange charges in 
the droplet. This has been done in the following way. 
From the solution of hydrodynamic equations 0-0 and 
(tT^ we know five quantities for each fluid: T™, 

and . In order to interpret them in thermody¬ 
namic terms, we should first determine six quantities for 
each fluid: ria, £a, Pa and 3 components of the 4-velocity 
Mq, proceeding from above five hydrodynamic quanti¬ 
ties, see Eqs. m and dJ. Naturally, five equations 
are not enough for determining six quantities. There¬ 
fore, we add one more equation to this set—the EoS. 
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Thus, the resulting thermodynamic quantities turn out 
to be EoS dependent. In particular, this scheme with 
the gas EoS results in a change of the hydrodynamic 4- 
velocity ('«fQ(gas)) as compared to that calculated with 
the EoS used in the hydrodynamic simulation. The sec¬ 
ond, already conventional step consists in determination 
of temperature and chemical potentials proceeding from 


baryonic, strange and energy densities and pressure. 

This is a kind of “freeze-out shock”, which however is 
completely different from that induced by the cut-Jiittner 
recipe. Now in terms of frozen-out droplets of various a- 
fluids, the spectrum of observable hadrons of a species 
with Cj baryon and e“ strange charges can be expressed 
as follows 


E 


dNg 

(Pp 




(proper) 


^'M'“ia(gas) 


(27r)3 


exp { 


„a,,*“(gas) 


„a,/“(gas 


^ / 7 "*«(gas)| 


± 1 


(33) 


where pa is degeneracy of the a particle, is the 

volume of the ia-droplet in its rest frame, and the sum 
runs over all frozen-out droplets of all fluids. If the a 
species is a baryon, the upper sign (-I-) should be taken 
in the denominator, if it is a meson, the lower sign (—). 

The described “fragmentation” method of freeze- 
out precisely conserves energy-momentum and various 
charges. However, it is also not free from problems. The 
microscopic justification of this method is still lacking. 
The problem of returning frozen-out particles into the hy¬ 
drodynamic phase still persists. Certainly, further search 
for a reliable freeze-out procedure is needed. A more con¬ 
sistent way of performing the freeze-out, the method of 
“continuous emission”, was proposed in Ref. ^3 ■ Avoid¬ 
ing sharply defined freeze-out hypersurface, this method 
considers a continuous emission of particles from a finite 
volume, governed by their mean free paths. Its predic¬ 
tions at the level of observables differ from those based on 
the Cooper-Erye recipe Unfortunately, this method 
is very difficult for the numerical implementation. 


III. SIMULATIONS OF NUCLEUS-NUCLEUS 
COLLISIONS 

The strategy of our simulations is as follows. At the 
first step, we try to reproduce the stopping power ob¬ 
served in proton rapidity distributions by means of fine 
fitting of friction forces of the model. In principle, fric¬ 
tion forces are EoS dependent both through scalar den¬ 
sities EJ and because of medium modifications of cross 
sections. In spite of this, our friction is just schemati¬ 
cally estimated proceeding from vacuum proton-proton 
cross sections and therefore does not comply with the 
EoS used. This is so even for the simplest hadronic EoS, 
used in the present paper. The list of uncertainties rele¬ 
vant to the friction forces has been discussed in Subsect, 
cm In view of this, the strategy of fine fitting of fric¬ 
tion forces is quite reasonable. Note that the success of 
such fit is not obvious in advance because by means of 
two functions ^h{s) and ^q{s) of a single variable (in the 


present case, only ^h{s)), see Eq. (t^ . we fit a function 
of three variables: rapidity, incident energy and impact 
parameter. It is worthwhile to mention that the fit of the 
friction is EoS dependent. It allows to reproduce proton 
rapidity distributions only with the particular EoS. In 
general, another EoS requires different fit. 

After the friction forces have been fixed, we tune the 
freeze-out energy density (in the reasonable range) in 
order to reproduce transverse-mass proton spectra and 
multiplicities of produced pions at comparatively low in¬ 
cident energies Ei^h < 20A GeV. At higher incident en¬ 
ergies, the pion multiplicity turns out to be weakly sen¬ 
sitive to variation of the freeze-out energy density. In 
principle, the freeze-out criterion could be different for 
different particle species and even for the chemical and 
thermal freeze-out. However, we keep it unique for all 
the cases in order to avoid multiplication of fitting pa¬ 
rameters. Therefore, fit of these quantities by means of 
a single parameter is not an obvious task. At higher in¬ 
cident energies Ujab > 30A GeV the pion multiplicities 
become sensitive to the formation time r. Therefore, the 
next step consists in tuning r (again in the reasonable 
range) to reproduce them already at higher incident en¬ 
ergies. These subsequent steps are simplified by the fact 
that proton rapidity distributions are only slightly sensi¬ 
tive to variation of r and the freeze-out energy density. 
After all these steps, all the model parameters got fixed, 
and all further calculations give pure predictions of the 
model. Our final aim is to find a EoS which reproduces 
in the best way the largest body of the observables of 
nuclear collisions in the incident energy range Ejab — (1- 
160)A GeV. 


A. Hadronic EoS 

We start our simulations with the purely hadronic EoS. 
This EoS, which was originally used in 2-fluid simulations 
[1(1 ITlL Il3 , was proposed in ^3 ■ It is a natural refer¬ 
ence point for any other more elaborate EoS. The energy 
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density and pressure are constructed as follows: 
e{nB,T) = egas{nB,T) + W{nB), (34) 

PinB,T) = P^^{nB,T)+nB^^^^^-W{nB),{^5) 

where £gas(?^s,ir) and PgasinB,T) are the energy den¬ 
sity and pressure of relativistic hadronic gas, respec¬ 
tively, which depend on baryon density ub and temper¬ 
ature T. The only difference from the ideal gas is that 
baryons are affected by a mean field [/(n^), i.e. the en¬ 
ergy of the a-baryon of mass Ma with momentum p is 
ea = (P^ + -I- baU{nB), where ba is the baryon 

number of the a-particle, and the potential U^ub) is 
parametrized as follows 


UijiB) = rriM 



-I- c(7 -t- 2) 



.(36) 


It depends only on the density ub- The self-consistent 
potential contribution to the energy density, IT (ns), is 

pTlB 

W{nB) = / U{n) dn. (37) 

Jo 

Parameters b, c and 7 are determined from the condi¬ 
tion that the cold nuclear matter saturates at no = 0.15 
fm“^ and e(no,T = 0)/no — mjv = —16 MeV, and in¬ 
compressibility of this nuclear matter is K = 210 MeV. 
Parametrization ll^ results in superluminal sound ve- 



FIG. 1: (Color online) Baryon-density dependence of the 
pressure at T = 0. Shaded region is the constraint derived 
from experimental data. Solid and dashed lines present the 
pressure P(ns,T = 0) and that for the ideal hadronic gas, 
Pgas(ns,T = 0) (cf. Eq. Idhll l. respectively. 


locity at densities ns/no > 8 . To preserve causality at 
high Ub, the following form of the energy density 


e{nB,T = 0) = nomjv 


A 



+ C + P 



is used at ubIuq > 6 . Parameters B and C are de¬ 
termined on the condition that e{nB,T = 0 ) and its two 


first derivatives are continuous at ub/tiq = 6 . As seen 
from Fig. Q] the pressure of this hadronic EoS is within 
the constraint given by Danielewicz et al. extracted from 
the analysis of flow of nuclear matter the AGS in¬ 

cident energies. 


B. Summary of Parameters: Friction, Freeze-Out, 

etc. 

For the sake of complete account, in this section we 
would like to summarize the parameters used in the 
present simulations. The specific reasons for choosing 
these parameters will be discussed in subsequent Subsec¬ 
tions. 

• The key quantity is the EoS, which was taken in the 
simple form of purely hadronic EoS, see Sect. ehi 
I n this EoS, 48 different hadronic species are taken 
into account. Each hadronic species includes all 
the relevant isospin states, e.g., the nucleon species 
includes proton and neutron. 

• The friction between baryon-rich fluids was fitted 
to reproduce the stopping power observed in pro¬ 
ton rapidity distributions, see Sect. irrm For this 
purpose the original friction, parametrized through 
experimental inclusive proton-proton cross sections 
| 22 j, was enhanced by means of the tuning factor 
ih 


ihi^) — 7h + 


In (s/( 2 mAr) 2 ) 


1 / 2 ]1/4 


(39) 


with 7 ^ = 1 and Ph = 0.75, see Fig. |21 The 
factor is not applicable here because of the pure 
hadronic nature of the EoS. 



FIG. 2: (Color online) Fitted friction enhancement for the 
pnrely hadronic EoS 14^ used in simulations, see Eq. lUB- 















The parameter of freeze-out energy density was fit¬ 
ted to reproduce transverse-mass proton spectra, 
see Subsect. ITTT PI and multiplicities of produced 
pions at comparatively low incident energies Tliab < 
20A GeV, see Subsect. IIII El It has been taken 
eirz — 0.2 GeV/fm^ at all incident energies with 
the only exception: Sbz ~ 0.1 GeV/fm^ at Euh < 
2A GeV. These e^z values are presented as approx¬ 
imate quantities because of the numerical realiza¬ 
tion of the freeze-out procedure (cf. App. A). 

The formation time of the fireball fluid was also fit¬ 
ted to reproduce the pion yield at higher incident 
energies ifiab > 30A GeV, see Subsect. 111111 It has 
been taken r = 2 fm/c. Note that r is the forma¬ 
tion time of a fluid element, consisting of a number 
of particles, which are generally not thermalized. 
This implies that the formation time of a separate 
particle is certainly shorter than r. 


Flab, A GeV 

1 

2 

4 

6 

8 

> 9 

^coal 

1.5 

1.2 

0.8 

0.4 

0.1 

0 


TABLE I: Coalescence parameters Gcoal used for simulations 
of Au+Au and Pb+Pb collisions at various incident energies 

.^lab- 

In order to compare with available experimental data, 
impact parameters used in calculations were taken either 
from experimental works, if they were evaluated there, 
or estimated based on the experimentally declared per¬ 
centage of the total reaction cross section, corresponding 
to a particular event selection. 


C. Proton Rapidity Distributions: Observable 
Stopping Power 


The contributions of strong decays of hadronic res¬ 
onances (i?) into spectra of stable hadrons is taken 
into account, as it is described in App. In the 
present calculation all spectral functions were taken 
without width 

Ar{s) =d{s-ml), (40) 

where rriR is the mass of the R resonance. 

For the reproduction of stable hadron multiplicities 
it is also important to take into account contribu¬ 
tions of weak decays with small values of cxu,, where 
Tu, is the inverse width of the decay. Feed back 
from weak decays of hadro ns wi th ct^ < few cm, 
i.e. A, E, A and E |l08| have been taken 

into account in yields of stable particles. In par¬ 
ticular, this is important for reproduction of the 
pion multiplicities (Subsects IIII El and IIII HU and 
the proper normalization of proton spectra (Sub¬ 
sects imn and IIII DI . On the other hand, weakly 
decaying hadrons with cr^, of the order of few me¬ 
ters, i.e. Ar+ and K~, were treated as stable 

particles. 

Light fragment formation (deutrons, tritons, ^He 
and '^He) is taken into account in terms of the co¬ 
alescence model, which is similar to that in Ap¬ 
pendix E of Ref. . In Ref. , the coalescence 
coefficients were fitted at the incident energy 0.8A 
GeV. To accommodate this formulation to differ¬ 
ent incident energies, we change only overall scale 
of the coalescence (the coefficient Ccoai) which ei¬ 
ther enhances (Ccoai > 1) or reduces (Ccoai < 1) its 
strength, keeping ratios of yields of various frag¬ 
ments the same as in M , see Table n At Ciab > 
9A GeV we do not apply any coalescence to calcu¬ 
lations of nucleon observables, since the respective 
correction is negligible. 


The nucleon rapidity distribution basically reflects the 
stopping power achieved in a nuclear collision. It is de¬ 
fined as 


dN]\[ 

dy 


cPpT 


E 


dNN ^ 

(pp 2^ (pp 


(41) 


in terms of frozen-out spectra of nucleons EdN^;/(Pp, 
see Eq. (ESI, and baryonic resonances. The second term 
under the integral takes into account the contribution of 
these resonance decays, cf. Eq. (IB5II . App. m the sum 
runs over all baryonic resonances R. Integration runs 
over the transverse (with respect to the beam) momen¬ 
tum pt ■ 

Either identified protons or difference between positive 
and negative hadrons are experimentally measured. The 
latter is associated with the proton-antiproton difference, 
{p — p). Since particles are not isotopically distinguished 
in our model, we estimate the proton distribution simply 
as 


dNp _ Z dNN , , 

dy A dy ' 

where Z and A are the proton and mass numbers, respec¬ 
tively, in the colliding (identical) nuclei. This estimate is 
quite reasonable at comparatively low incident energies. 
At higher energies, when abundant particle production 
starts, this recipe somewhat underestimates the proton 
number. The reason is that newly produced particles 
tend to restore the isotopic symmetry in the system and 
hence increase the number of protons as compared with 
its initial value. The (p—p) quantity is much less sensitive 
to the effect of newly produced particles, since their con¬ 
tribution is essentially canceled in the difference (p — p). 
Therefore, the Z/A scaling of the nucleon-antinucleon 
difference is a reasonable approximation for the (p — p) 
quantity even at high incident energies. As mentioned 
above, our aim is to reproduce these distributions in a 























dN/dy dN/dy 
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wide range of incident energies from lA to about 160A 
GeV. 

In general, nucleon observables are more robust to vari¬ 
ations of physical parameters than other probes, since 
they are essentially confined by the baryon number con¬ 
servation. Our first observation is that the nucleon 
and, in particular, proton rapidity distributions are only 
weakly sensitive to variation of the freeze-out energy den¬ 
sity Efrz and the formation time of the fireball fluid r. As 
for Efi-z, there is a certain compensation between effects 
of collective motion and internal excitation. Both these 
effects produce similar consequences in nucleon spectra. 
For instance, if we allow the system to evolve longer (the 



-1 -0.5 0 0.5 1 

y/Ycm. 


FIG. 3: (Color online) Rapidity spectra of protons (upper 
panel) and (p — p) (lower panel) from central heavy-ion col¬ 
lisions. Bold solid lines and bold dashed line (upper panel) 
correspond to 3-fluid hydrodynamic calculations at different 
impact parameters. Thin dashed-dotted and short-dashed 
lines are the appropriate UrQMD and HSD model results Id^ . 
respectively. Full symbols display measured experimental 
points, whereas the open ones are those reflected with respect 
to the mid rapidit y p oint. Experimental data are taken by col¬ 
laborations E802 lii, E877 E917 [H, NA49 [5^.|54|5^. 

The percentage shows the fraction of the total reaction cross 
section, corresponding to experimental selection of central 
events. 


lower Efrz), the collective motion of the matter becomes 
more developed while its internal excitation drops down 
(i.e. the matter cools down). This counteraction results 
only in slight changes in nucleon rapidity distributions. 
As for the formation time, the density of the produced 
fireball fluid is not high enough even at the highest con¬ 
sidered energy. Therefore, its interaction with bar yon - 
rich fluids only slightly affects the baryon subsystem . 
This weak dependence on ef^z and r allows us to use nu¬ 
cleon rapidity distributions to fit the observable stopping 
power, to which nucleon observables are indeed sensitive. 

The second main conclusion is that the original fric¬ 
tion between baryon-rich fluids, estimated in Ref. 
proceeding from free proton-proton cross sections, is ev¬ 
idently insufhcient to reproduce observable stopping of 
nuclear matter. Therefore, this friction is htted to re¬ 
produce this stopping power observed in proton rapidity 
distributions. For this purpose the original friction was 
enhanced by means of the tuning factor cf. Eq. inm . 
Thus, the friction enhancement is the larger, the higher 
incident energy is. At the lowest considered energy of lA 
GeV, the friction turns out to be approximately 2 times 
enhanced. This is similar to earlier results of the 2-fluid 
model with mean mesonic helds M- There it was found 
out that for the proper reproduction of data on heavy-ion 
collisions in the energy range from 0.4 to 0.8A GeV the 
enhancement factor of 3 was required. 

Fig. 0 illustrates the overall quality of reproduction 
of experimental data by our hydrodynamic calculations 
as well as by other transport simulations for two basic 
incident energies of 10 and 160A GeV. The 3-fluid simu¬ 
lations were performed at impact parameters which cor¬ 
respond to the experimentally declared fraction of the 
total reaction cross section, accumulated in the experi¬ 
mentally selected central events, assuming sharp cut off 
in impact parameters. The side bumps of the HSD calcu¬ 
lations correspond to spectator parts of colliding nuclei. 
In our model, these spectator parts were cut off, based 
on simple criterion: the energy per baryon is less than 
the nucleon mass, which is the case, when a nucleon is at 
least loosely bound in the matter. As seen, experimen¬ 
tal data from different experiments somewhat differ, and 
therefore one should not expect their reproduction in the 
model calculations within better than ~ 10%. Uncer¬ 
tainties in the impact parameter at the level of about 0.5 
fm do not noticeably affect hydrodynamic rapidity spec¬ 
tra. In general, there is a reasonable agreement between 
hydrodynamic and kinetic calculations. 

Comparison with experimental data for identified pro¬ 
tons is presented in Fig. 0for the SIS-AGS energy range, 
where various selections of non-central nuclear collisions 
are also considered. Here and below, impact parame¬ 
ters for each set of non-central interactions were chosen 
accordingly the experimentally declared fractions of the 
total reaction cross section, accumulated in the set of ex¬ 
perimentally selected events, assuming sharp cut off in 
impact parameters. 

At the beam energies E\ah of the order of lA GeV, a 
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FIG. 4: (Color online) Proton rapidity spectra at SIS and AGS energies for various impact parameter. Solid lines correspond 
to 3-fluid calculations. The dashed line is that for deuteron production. For comparison the 1-fluid result is shown by the 
dotted-dashed line at i?iab = 1.06A GeV. Experimental points at different energies are taken from at Ejab = 1.06A GeV, 
at 2A and 4^4 GeV, and at 6A, 8A and 10.5A GeV. The percentage indicates the fraction of the total reaction cross 
section, corresponding to experimentally selected events. 


large fraction of nucleons is produced with low relative 
velocities. Therefore, they may coalesce forming light 
fragments. The coalescence formulation accepted here is 
similar to that in Appendix E of Ref. M Only light 
fragments (deutrons, tritons, ^He and '^He) are taken into 
account. To fit this formulation to different incident ener¬ 
gies, we change only overall scale of the coalescence (the 
coefficient Ccoai) which either enhances (Ccoai > 1) or 
reduces (Ccoai < 1) its strength, keeping ratios of yields 
of various fragments the same as in [^, see Tabled 

At the SIS energy we describe reasonably well both 


proton and deuteron rapidity spectra. The fact that the 
experimental distribution is somewhat wider than the 
calculated one can be explained by that the set of experi¬ 
mentally selected events in fact contains a certain admix¬ 
ture of semi-central and peripheral events, whereas the 
calculation was performed for the single central impact 
parameter. Indeed, at comparatively low energies the re¬ 
liable selection of central events is problematic . This 
non-central admixture makes the distributions wider as 
compared to what it would be for the perfect central se¬ 
lection. At this energy the stopping power is rather high. 
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FIG. 5: (Color online) (p — p) rapidity distribution from cen¬ 
tral Pb-|-Pb collisions at SPS energies. Solid lines are 3-fiuid 
calculations, whereas dashed-dotted and dashed lines are the 
corresponding predictions of the HSD and UrQMD models 
m, respectively. NA49 experimental data are from Ref. 
for i^iab = 40A and 80A GeV and from Refs. 1^ 1^ for 
^lab = 158A GeV. The percentage indicates the fraction of 
the total reaction cross section, corresponding to experimen¬ 
tally selected events. 


This is seen from both the Gaussian-like shape of the ra¬ 
pidity distributions and also from the fact that the 3-fluid 
results are quite close to those of the conventional 1-fluid 
calculations (see Fig. El). A separate code was used for 
these 1-fluid calculations. With the increase of the inci¬ 
dent energy, the spectrum shape starts to differ from the 
Gaussian one, getting more and more flat at the mid ra¬ 
pidity. The 3-fluid model reasonably reproduces the de¬ 
pendence of the spectra on both the incident energy i?iab 


and the impact parameter b. Note that no extra tuning 
of normalization of the spectra was done. The values of 
impact parameters were taken from the experimental es¬ 
timate of centrality of nuclear interactions . Though 
there are some uncertainties in this estimate, we applied 
no special tuning to values of these impact parameters. 

The (jj — p) rapidity spectra at the SPS energies are 
shown in Fig. El A minimum of hydrodynamic distribu¬ 
tions dN/dy in the mid rapidity region complies with ex¬ 
perimental observations at Euh = 158A GeV. It survives 
at lower incident energies, up to Aiab = 40A GeV. This 
is in contrast to simulations based on kinetic transport 
codes 0 , which predict flat or even peacked distribu¬ 
tions in the mid rapidity region at £^iab = 40A and 80A 
GeV. Available experimental points here do not allow us 
to verify these different predictions. Note that experi¬ 
mental points at Ei^h = 40A and 80A GeV correspond 
to identified protons, while all presented calculations are 
related to {p — p). 


D. Proton Transverse Mass Distributions 

Once the friction has been already fitted to reproduce 
proton rapidity distributions, we can vary only the freeze- 
out energy density e^z and the formation time of the 
fireball fluid r. In contrast to rapidity distributions, the 
proton transverse-mass distributions turn out to be more 
sensitive to since their slopes reflect the effective 
temperature of the frozen-out distributions. This is be¬ 
cause the internal excitation certainly dominates in the 
transverse direction as compared to the collective motion. 
At the same time, these distributions are rather insensi¬ 
tive to r, as well as all other baryonic (however, not an- 
tibaryonic) quantities. Therefore, to reproduce them, we 
choose Efrz — 0.2 GeV/fm^ at all incident energies with 
the only exception: er^z ~ 0.1 GeV/fm^ at Ei^h < 2A 
GeV. As for low energies, Ei^h < 2A GeV, we had to 
reduce euz, since otherwise the freeze-out would occur at 
the very early stage of the collision. The approximate 
nature of the e^z values results from the numerical real¬ 
ization of the freeze-out procedure (cf. App. A). 

Proton transverse-mass spectra at the AGS energies 
are exemplihed in Fig. El They exhibit a typical expo¬ 
nential fall-off. As seen, the 3-fluid model well reproduces 
both the normalization and the slopes of this fall-off, as 
well as their rapidity dependence. The right panel of Fig. 
El shows the dependence of the transverse mass distribu¬ 
tions at the mid rapidity on the impact parameter b. As 
before, the values of impact parameters were taken from 
the experimental estimate of centrality of nuclear inter¬ 
action M- The impact-parameter dependence is reason¬ 
ably reproduced by the model as well even at relatively 
large impact parameters. 

The beam-energy dependence of the proton transverse- 
mass distributions is presented in Fig. [3 for central 
Pb-I-Pb collisions at the SPS energies. Overall, these 
distributions are in agreement with experimental data. 
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FIG. 6: (Color online) Transverse mass distributions of protons from Au(8A GeV)+Au collisions at & = 2 fm and different 
rapidities in the c.m. system (left panel) and at the mid rapidity and different impact parameters (right panel), cf. the right 
panel of Fig. 0] For clarity of representation, every next data set and the corresponding curve (from top to bottom) is multiplied 
by the additional factor 0.1. Experimental data are from m. 



mT-m [GeV] 

FIG. 7: (Color online) Transverse mass distributions of pro¬ 
tons at the mid rapidity from central Pb-|-Pb collisions at 
incident energies = 158A, 80A and 40A GeV and im¬ 
pact parameter b — 2.5 fm, cf. Fig. El NA49 experimental 
data are taken from (squares), (open circles) and l5^ 
(full circles). 



y 

FIG. 8: (Color online) Pion rapidity spectra from central 
Au(10.5A GeV)+Au collisions. Bold solid and dashed lines 
represent 3-fluid calculations at fc = 2.5 and 2.0 fm, respec¬ 
tively. Thin dashed-dotted and short-dashed lines correspond 
spectra from kinetic simulations within HSD and UrQMD 
models ll^l . The •n~ (full squares and triangles) and vr'*' (full 
circles) data are measured by E895 and E877 IHlt collabo¬ 
rations. The percentage indicates the fraction of the total re¬ 
action cross section, corresponding to experimentally selected 
events. Open symbols are obtained by reflecting the full ones 
with respect to the mid rapidity. 


However, the 3-fluid model does not exhibit deviation 
from the exponential fall-off at (mp — m) <0.2 GeV, ob¬ 
served in experiment, most clearly at = 158A GeV. 
The same feature of the hydrodynamic calculation was 


earlier reported in Ref. [^. As it was shown [^, a 
post-hydro kinetic evolution is required to produce the 
observable two-slope structure of the WT-spectra. In 
our model such a post-hydro evolution is absent. It is 
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FIG. 9: (Color online) Pion rapidity spectra from central collisions at the AGS (left panel) and SPS (right panel) energies. 
Experimental points are taken from Im (AGS energies) and IH^I (SPS energies). Notation is the same as in Fig. 


worthwhile to note that the above problem is not an in¬ 
alienable feature of any hydrodynamic calculation. For 
instance, the deviation from the exponential fall-off was 
reproduced in calculations by Kolb et al. . 


E. Pion Rapidity Distributions 

Production of new particles is closely related to the 
amount of entropy accumulated in the system. At the 
late stage of the collision, the system expansion is isoen- 
tropic, i.e. the total entropy is conserved. This is indeed 
so in our simulations, as we have checked it. This fact 
implies that at high incident energies Eiah > 40A GeV 
the total number of produced pions is approximately con¬ 
served during this expansion stage, since thermal pions 
are the dominant component among the produced parti¬ 
cles. From the practical point of view, it means that the 


pion number is approximately independent of the freeze- 
out energy density efrz. This is indeed observed in actual 
simulations. However, at high incident energies the pion 
number depends on the formation time of fireball fluid 
T, since 20-^30% of the pions are produced in the fire¬ 
ball fluid. The shorter r is, the earlier the fireball fluid 
starts to interact with baryonic subsystem, and hence 
the fewer pions survive in this fluid. This implies that 
the formation time effectively tunes the strength of the 
fireball-projectile(target) friction which was only roughly 
estimated (cf. subsect. liTTl . 

At the first glance, the above speculation contradicts 
to the pion absorption mechanism (cf. in Subsect. irm . 
which assumes that pions are simply captured by the 
baryon-rich fluids without loosing their number. How¬ 
ever, the pion number is not a conserved quantity. In 
the baryon-rich fluids the energy of the captured pions is 
thermally redistributed (accordingly to thermal chemical 
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FIG. 10: (Color online) Rapidity spectra of A + E° hyperons (left panel) and A + E° antihyperons (right panel) from central {b 
= 2.5 fm) Pb+Pb collisions. Solid lines represent results of the 3-fluid model. Contributions from the hreball fluid are shown 
by dashed lines. Preliminary experimental data are taken from Ref. |^. The percentage indicates the fraction of the total 
reaction cross section, corresponding to experimentally selected events. 


equilibrium) between kinetic energy of (mainly) baryons 
and thermally produced (mainly) pions. Therefore, their 
number is effectively reduced, since in the baryon-free 
fireball fluid the same energy was mainly accumulated in 
thermal pions. 

At lower incident energies Aiab < 20A GeV, the ther¬ 
mal pion production is not already so dominant but pi¬ 
ons are rather produced through decays of resonances. 
Therefore, their number starts to depend on ef^z- The 
later freeze-out occurs, the fewer pions are produced. 
At the same time, the contribution of the produced 
fireball fluid becomes less important, and hence the t- 
dependence of the pion multiplicity becomes weak. At 
intermediate energies 20 < Aiab < 40A GeV, there is 
a moderate dependence of the pion multiplicity on both 
Efrz and r. Note that these parameters, e^z and r, mostly 
affect the overall normalization of pion rapidity spectra. 


whereas the shape of these spectra mainly depend on the 
stopping power, which is kept fixed here. 

As both the friction and £frz have been already fixed 
above, the pion rapidity distributions at low incident en¬ 
ergies Aiab < 20A GeV can be considered as predictions 
of the model. In order to reproduce the pion yield at high 
incident energies ifiab > 30A GeV , the formation time 
of the fireball fluid has been taken r = 2 fm/c. Note that 
T is the formation time of a fluid element, consisting of 
a number of particles (cf. Eq. 11411 and discussion above 
it). This implies that the formation time of a separate 
particle (pion) at rest is certainly shorter than r. As¬ 
suming that the temperature scale of the fireball fluid is 
approximately 100 MeV, we estimate that the formation 
time of a separate pion at rest approximately equals 1 
fm/c. This is in good agreement with the estimate of 
the string-formation time. However, one should keep in 
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FIG. 11: (Color online) Antiproton rapidity spectrum from 
central Pb(158A GeV)+Pb collisions. Notation is the same 
as in Fig. ITUl Experimental data are taken from Ref. 1^ . 

mind that the latter value refers to a separate string in 
vacuum, i.e. without any effects of string interactions 
(junctions and color-rope formation). At the same time, 
our r has an in-medium sense with all these multi-particle 
effects included. Therefore, their coincidence should be 
taken with caution. Note that with r = 0 the multiplic¬ 
ity of pions (as well as other newly produced particles) 
would be 10-15% underestimated at high incident ener¬ 
gies (> 30A GeV). All other observables would remain 
approximately the same. 

Since particles are not isotopically distinguished in our 
model, we assume that numbers of pions of each charge 
are equal, i.e. = N^- = N^o = where is 

the calculated total number of pions. Of course, this is 
a rough estimate of 7r+ and tt~ yields, since always 
exceeds At,r+ because of the initial isotopic asymmetry 
of colliding nuclei. Therefore, only if our calculation of 
Nt^/Z complies with experimental -\- A^^-)/2, we 

refer to this as ” a good agreement”. 

Fig. IHlrepresents pion rapidity distributions in central 
Au-I-Au collisions at the incident energy E\ah = 10.5A 
GeV. Sensitivity of this distribution to the variation of 
the impact parameter is also demonstrated. As seen, the 
hydrodynamic results indeed fall in between the tt’’' and 
7r“ data, hence they comply with these data. In contrast, 
the Tr’*' rapidity spectra, calculated in two transport mod¬ 
els, HSD and UrQMD models 0 also displayed in Fig. 
IHl closely follow experimental points but for negative pi¬ 
ons. The overestimate of the pion yield by the kinetic 
models can be attributed to a lack of collective interac¬ 
tion in the EoS’s corresponding to those models. 

All above said is in fact true for the whole energy range 
under consideration, as it is demonstrated in Fig. 

The 3-fluid model reasonably reproduces the pion dis¬ 
tributions, while the HSD and UrQMD models certainly 
overestimate them. 


F. Rare Particle Production 

Now all the parameters of the model are fixed. There¬ 
fore, all further calculations can be treated as predictions 
of the 3-fluid model. 

As an example of the hydrodynamic description of rare 
channels, the rapidity spectra of A and A hyperons are 
presented in Fig. ^1 Contributions from decays of higher 
resonances are taken into account. Both shape and abso¬ 
lute value of hyperon and antihyperon spectra are repro¬ 
duced surprisingly well. The A hyperons originate mainly 
from baryon-rich fluids and shapes of their spectra closely 
follow those of the protons. Although the experimental 
data are preliminary, one may note that agreement with 
experiment becomes certainly worse at the incident en¬ 
ergy Aiab = 40A GeV. Unfortunately, the only measured 
point in proton rapidity distribution (see Fig. O does 
not allow us to conclude on either similarity or difference 
of the A and p distributions at this energy. In contrast, A 
antihyperons are dominantly created in the baryon-free 
fireball and have typical single-bump thermal spectra. 

At the same time the experimental data on antipro¬ 
ton production are certainly overestimated in the our 
model, as illustrated in Fig. Since particles are 

not isotopically distinguished in our model, we have esti¬ 
mated the antiproton yield as half of that of antinucleon: 
Np = assuming that pp and nn pairs are produced 

with approximately equal probability. The contributions 
of weak-decay channels A ^ N + tt and E —> V -|- tt are 
excluded. 
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FIG. 12: (Color online) Transverse-momentum flow of nu¬ 
cleons and protons as a function of rapidity for semi-central 
Au+Au collisions at — 10.5A GeV. The 3-fluid calcu¬ 
lations are presented for 6 = 4 fm (solid line) and 6 = 6 fm 
(dashed line). The data for protons (circles) and all baryons 
(triangles) are taken from Ref. Full symbols correspond 
to measured data, open ones are those reflected with respect 
to the mid rapidity point. 
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G. Flow 

Flow quantities of different types quantify space- 
momentum correlations of collective motion of the 
strongly interacting matter. This collective motion is es¬ 
sentially caused by the pressure gradients arising during 
the evolution of the collision and hence is intimately re¬ 
lated to the EoS and, in particular, to a possible phase 
transitions. A spectacular loss of correlation between the 


observed particle transverse momenta and the reaction 
plane, which gives rise to dramatic reduction of the di¬ 
rected flow, has been predicted by Rischke et al. M 
in the conventional hydrodynamic model with the bag- 
model EoS. The subsequent studies showed that the ob¬ 
served signal essentially depends not only on the EoS but 
also on the collision dynamics [Ulli. 

The conventional transverse-momentum flow of an a 
species is defined as cf. Eq. gu, 




J tfpT Px 

E dNa/d^p + '^E /d^p 

R 

J d'^PT 

E dNa/d^p + ^ A 

R 



(43) 


where px is the transverse momentum of a particle in the 
reaction plane, and integration runs over the transverse 
momentum px- The second term in square brackets takes 
into account the contribution of resonance decays, result¬ 
ing in a production, cf. App. m the sum runs over all 
relevant resonances R. 

As seen in Fig. im the 3-fluid model reasonably repro¬ 
duces a general trend of the baryonic {px)(jj) distribution 
at Aiab = 10.5A GeV, exhibiting two clear peaks near 
the target and projectile rapidities. In fact, we compute 
{Px)[y) of so-called primordial nucleons, which later may 


coalesce, forming light fragments. In view of this, it is 
not surprising that agreement with the nucleon {px){y) 
data, which indeed take into account contribution of light 
nuclear fragments, seems to be better than that with 
identilied-proton data. 

At high energies the azimuthal asymmetry is usu¬ 
ally characterized by the first and second coefficients of 
Fourier expansion of the azimuthal-angle dependence of 
the single-particle distribution function, i.e. by the di¬ 
rected flow vi = (cos <j)) and elliptic flow V 2 = (cos 2(j)): 


v^i\y) 


4“^ (y) 


J d^priPx/pT) 

E dNa/d^p + ^ A /d^p 

R 

J 

d^PT 

E dNa/d^p + '^E /d^p 

R 



J d^PT [{pI-pI)/p^ 

E dNa/d^p + '^E /d^p 

R 

J d^PT 

E dNa/d^p + '^E /d^p 

R 



(44) 


(45) 


Examples of the vi and V 2 flow for protons and pions at 
Aiab = 40A and 158A GeV are presented in Figs and 
d The declared experimental percentage of the total 
reaction cross section, related to the experimental event 
selection, allows us to estimate the corresponding impact 
parameter as b = 5.6 fm. In view of uncertainty of this 
estimate and in order to reveal the model dependence on 
the impact parameter, we also show calculations with b = 


4 fm. Flows of the fireball fluid calculated at 6 = 4 and 
5.6 fm are hardly distinguishable. Several sets of data, 
which noticeably differ from each other, are shown in 
these figures. The ’’standard” method for evalu¬ 

ating the flow coefficients requires an event-by-event esti¬ 
mate of the reaction plane with which outgoing hadrons 
correlate. However, this method does not discriminate 
other sources of correlations, like those due to global mo- 
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FIG. 13: (Color online) Directed flows of protons (upper panels) and charged pions (lower panels) in semi-central Pb-|-Pb 
collisions at Slab = 158^ (right panels) and 40A (left panels) GeV as functions of rapidity. The 3-fluid calculations are 
presented for b — 5.6 fm (solid line) and b = 4 fm (dashed line). The dot-dashed lines show contributions of the fireball fluid. 
Experimental data are taken from Ref. 1^ . These data were obtained by two different experimental procedures: the standard 
one (ni, 2 (st)) and the method of n-particle correlations (vi^ 2 (n)). Full symbols correspond to measured data, while the open 
symbols are those reflected with respect to the mid rapidity. Updated data of the NA49 Collaboration (ui, 2 (st) — 98) with 
acceptance 0.05 < pr < 0.35 GeV/c for pions and 0.6 < pr < 2.0 Gev/c for protons, are also shown. 


mentum conservation, resonance decays, etc. Recently a 
new method of n-particle correlations has been proposed, 
which allows to get rid of these non-flow correlations in 
extracting vi and V 2 from genuine azimuthal correlations 
1^ . Note that only statistical errors are indicated in 
these hgures. The systematic error for protons is 0.005 
for V 2 and 0.01 for mat 158A GeV, they can be by 50% 
larger at 40A GeV |^ . 

At these energies the 3-fluid vi flow significantly differs 
from the data. Overall, the 3-fluid model predicts an es¬ 
sentially stronger directed flow than that experimentally 
observed. This disagreement can not be caused by nei¬ 
ther uncertainties in the impact parameters nor applica¬ 
tion of different methods for measuring the directed flow 
1^ . Moreover, the calculated pion directed flow closely 
follows the pattern of the proton one, while the pion vi 
data reveal anticorrelation with proton flow. Note that 
the vi flow in the fireball fluid is very weak. Therefore, it 
is not surprising that pion and proton vi are correlated 
in the 3-fluid model, since they reflect the same collective 
motion of the baryon-rich fluids. 

The reasons for this poor reproduction of vi can be 
two-fold. Certainly the first reason consists in disre¬ 


garding the fact that a part of frozen-out particles is 
“shadowed” by still hydrodynamically evolving matter 
(cf. discussion in Subsect. IllUll . This shadowing means 
that frozen-out particles cannot freely propagate through 
the region still occupied by the hydrodynamically evolved 
matte r bu t rather get reabsorbed into the hydrodynamic 
phase [T?)^ . Apparently, the baryon directed flow is less 
affected by this shadowing. The reason is that the baryon 
directed flow reveals the collective flow of matter, which is 
mainly built of baryons as the most abundant and heavy 
component of the system. This collective flow is mainly 
formed at the early stage of the reaction. Baryon rescat¬ 
terings within this earlier-formed collective flow at later 
(freeze-out) stages do not essentially alter the collective 
transverse momentum of the matter. At the same time, 
the pions can be much stronger affected by this shadow¬ 
ing, since they are screened by the predominantly bary- 
onic matter, where pions may be essentially decelerated 
or even absorbed. This can drastically change the pion 

Vi. 

The second possible reason is very probable from our 
point of view. It could happen that the hadron EoS 
used in these simulations is too hard, i.e. results in too 
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FIG. 14: (Color online) The same as in Fig. 1 141 but for the elliptic flow. 


strong bounce-off of matter. A softer Eo^ in particu¬ 
lar due to possible phase transition , would re¬ 

duce the strength of the proton directed flow in favor 
of its better reproduction. Then the problem of corre¬ 
lation/anticorrelation with the pion flow would be more 
delicate, because the shadowing effect could easily change 
the sign of the pion flow. An additional argument in favor 
of a softer EoS is the success of microscopic RQMD 
and UrQMD models, which correspond to essentially 
softer EoS, since it is closer to that of the gas. These mod¬ 
els qualitatively reproduced early measurements of vi (y) 
and V 2 {y) for both protons and pions at Eiab = 158A 
GeV. Recently a good microscopic description of the dif¬ 
ferential flow at 40A GeV was also obtained [7n|. 

The calculated elliptic flow of protons (see Fig. nui is 
positive (”in plane”) and somewhat overestimates exper¬ 
imental points at both incident energies, 40A and 158A 
GeV. Though data at 40A GeV are still controversial, this 
overestimate again indicates that the used hadronic EoS 
is too hard. At the same time, the hydrodynamic pion 
elliptic flow is in surprisingly good agreement with ex¬ 
perimental data of the NA49 Gollaboration both in mag¬ 
nitude and shape, exhibiting a shallow minimum at the 
mid rapidity. The elliptic flow V 2 results from the initial 
spatial asymmetry of non-central nucleus-nucleus colli¬ 
sions. The overlap lens-shaped geometry of two nearly 
thermalized nuclei is reproduced correctly in the 3-fluid 
model. As this lens-shaped matter expands, it produces 
the elliptic flow. Note that mesons emitted from the fire¬ 


ball fluid have negative (’’out of plane”) V 2 flow. 

The earlier pion V 2 data of the NA49 collaboration, 
taken at smaller acceptance at the top SPS energy [^ . 
have been already analyzed within hydrodynamic ap¬ 
proaches. In the expansion model the Bjorken scaling 
solution 1^ was assumed for longitudinal evolution and 
2D hydro was solved numerically for transverse one. In 
this way, the elliptic flow could be estimated only at the 
mid rapidity point. The full 3D expansion model with 
postulated initial conditions was applied to the meson 
elliptic flow by Hirano |7^ . In a qualitative agreement 
with Ref. El , it was found that p-meson decays result in 
almost vanishing azimuthal anisotropy of pions near the 
mid rapidity. Let us remind that decays of all relevant 
resonances are taken into account in our model. 

Recently, collective flows in heavy ion collisions from 
AGS to SPS energies were systematically studied in a 
transport model with various assumptions on the nuclear 
mean field M- It was found that momentum depen¬ 
dence in the nuclear mean field is of prime importance 
for the reproduction of the vi and V 2 flows. It turned out 
that generally our results are rather close to those of Ref. 
|74| . if the mean-field momentum dependence is absent. 
This is precisely the case in our EoS. It is of interest to 
note that similarly to our results but in contrast with 
the experiment, the pion directed flow at Eiab = 158A 
GeV correlates with the proton v-\ independently of the 
version chosen for the mean field l74l. 
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FIG. 15: (Color online) Multiplicities for central Au+Au (at AGS energies) and Pb+Pb (at SPS energies) collisions as functions 
of the incident energy. O pen triangles connecte d by solid lines represent the 3-fluid results. The experimental data (full circles) 
are taken from Refs. IH^. l75l. It^ IttI 173. It^ l8ll 1^ 1^ 1^ 1^ for SPS energies and from Refs. for AGS 

energies. 




FIG. 16; (Color online) The same as in Fig. ^1 but for multiplicities of multi-strange hyperons. The experimental data are 
taken from Refs. 1^ 194.19^ : full circles for hyperons, and full squares for antihyperons. 


H. Multiplicities 

Though the quality of reproduction of particle multi¬ 
plicities was already clear from the above presented ra¬ 
pidity distributions, we would like to summarize them 
in this subsection, as well as present those for hadrons 
not mentioned above. In Fig. ESI particle multiplicities 
for central Au-|-Au (at AGS energies) and Pb-|-Pb (at 
SPS energies) collisions as functions of the incident en¬ 
ergy are presented and confronted with the experimen¬ 


tal data. The experimental data are taken from Refs. 

for SPS ener¬ 
gies and from Refs. for AGS energies. 

These data are presented as they were summarized in 
Ref. [qJ , where they were scaled to the same 5% experi¬ 
mental trigger. Because of this reduction to the same 5% 
trigger, all the simulations were performed at the same 
central impact parameters corresponding to this trigger: 
6 = 2 fm for Au-|-Au and b = 2.5 fm for Pb-|-Pb. Gor- 
rections for feed down from weak decays (mostly 
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and their antiparticles) were not applied in the original 
published data. They were, however, estimated to be 6% 
for A and 12% for A [t^. Based on this estimate the A 
and A yields given in Ref. were reduced by 6% and 
12%, respectively. This also complies with conditions of 
our simulations. 


We compare our calculated pion multiplicities with the 
experimental results for half-sum of those for 7r+ and 
7r“, since the model does not distinguish the isospin of 
particles. Note that the pion multiplicities were fitted 
by means of the freeze-out parameter and the formation 
time of the fireball fluid, while all other multiplicities are 
already predictions of the model. As seen, the model 
reasonably reproduces the energy dependence of various 
multiplicities, except that for K~. The latter is strange, 
especially in view of that even the rare-channel parti¬ 
cles, such as (j) and A, are reasonably reproduced without 
any additional tuning. Probably, a kind of post-h ydr o 
kinetic evolution, similar to that performed in Ref. |22l| . 
is required for a proper reproduction of the K~ multi¬ 
plicity. Indeed, cross sections of the reactions KN —> ttA 
and KN —> ttS are very high at low relative momenta 
The first reaction transforms a part of K~ into A. 
The resulting A enrichment would not contradict to data, 
since the NA57 data for the A production are cer¬ 
tainly above the NA49 data presented in Fig. The 
second reaction results in a loss of strangeness because of 
the weak decay ^ Ntt. In fact, this is the dominant 
channel for the K absorption on nucleons at low ener¬ 
gies. In our model such a post-hydro evolution is absent. 
It is worthwhile to note that the statistical-model fit of 
the data also overestimate K~ multiplicities at the SPS 
energies |^ . 


At comparatively low incident energies i?iab < lOA 
GeV, the grand canonical treatment of the strangeness 
production, used in the 3-fluid model, becomes poorly 
applicable. Therefore, we do not present the 3-fluid pre¬ 
dictions at these energies. At energies E\ah = (30-40)A 
GeV the observed multiplicities of and A are cer¬ 
tainly above the smooth curve predicted by the 3-fluid 
model. This correlates with already above observed fact 
(see Figs ^1 and d that agreement with data in this 
energy region is certainly worse than at other energies. 


The excitation functions for production of multi- 
strange hyperons, S and Sd, and the corresponding an- 
tihyperons are presented in Fig. d In Ref. 113, fl” 
and 0+ were not separated at i?iab = 40A GeV. Due to 
that their total yield is plotted as since the yield 
of D'’” antihyperons is much lower than that of It is 
seen that calculated multiplicities of multi-strange hyper¬ 
ons are also in a reasonable agreement with the available 
data. Note that these yiel ds are strongly underpredicted 
in the UrQMD model ' jsnll . 




FIG. 17: (Color online) Temporal evolution of average baryon 
(upper panel) and energy (lower panel) densities for central 
Pb-|-Pb collisions within 3-fluid model. The time is counted 
in the c.m. frame of the colliding nuclei. 


IV. GLOBAL EVOLUTION OF NUCLEAR 
COLLISIONS 


In the preceding Section, it was shown that the 3-fluid 
hydrodynamic model with the simplest, purely hadronic 
EoS is able to reasonably reproduce a great body of ex¬ 
perimental data on relativistic heavy-ion collisions. The 
considered observables characterize a state of the collid¬ 
ing system at the freeze-out stage. With thus fixed pa¬ 
rameters of the model, it is possible to address the ques¬ 
tion of global evolution of nuclear collisions, e.g., the val¬ 
ues of baryon and energy densities achieved in the course 
of them. An important question also is how long and in 
which volume these achieved values survive. Dynamics 
of heavy ion collisions at different bombarding energies 
is illustrated in Fig. El The shown baryon and energy 
densities are averaged over the whole space occupied by 
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FIG. 18: Time evolution of the entropy per baryon number 
in Pb+Pb central (6 = 2.5 fm) collisions at various incident 
energies £^iab- The total entropy is displayed by solid lines, 
while the entropy of baryon-rich subsystem - by dashed lines. 


fluids 


(ns) 

{£) 


dPx ub W[x)/ J dPx W{x) 
d^x eW{x)/ J d^x W(x) , 


(46) 

(47) 


where e and ub are the proper (i.e. in the local rest 
frame) densities, and the weight function is taken equal 
to the local proper baryon density W{x) = nB{x). This 
weight function was chosen because we are primarily in¬ 
terested here in evolution of dense baryonic matter. The 
3-fluid nonequilibrium is quite strong at the initial stage 
of the collision. The mean values of e and n b are calcu¬ 
lated by means of averaging these quantities correspond¬ 
ing to either separate fluids (if their mutual stopping has 
not occurred) or the unified fluid (if the full stopping 
has happened). To characterize the degree of stopping, 
we used the auxiliary function introduced in Eq. (P)l. 
which reveals quite a sharp transition between trans¬ 
parency (d ~ 1) and full stopping (d ~ 0). This conven¬ 
tion has been chosen in order to map the nonequilibrium 
configuration in terms of equilibrium quantities, avoid¬ 
ing unphysical contributions of the Lorentz contraction 
and collective motion into quantities (e) and (ns). How¬ 
ever, these results at early stages of collisions should be 
taken with some care because of certain ambiguity of this 
mapping convention. 

As seen from Fig. El the average baryon density 
(ns) increases quickly reaching a maximum, and then 
the expansion stage comes. The maximal compression 
of baryonic matter correlates with maximal overlap of 
colliding nuclei. Evolution of the average energy den¬ 
sity (e) proceeds in a similar way. Maximal values of 
baryon and energy densities are high enough, exceeding 
hypothetic threshold for deconfinement phase transition 


e ~ 1 GeVj jrn? even at ~ lOA GeV. Howeve^these 
values are slightly lower than those reported in Q], be¬ 
cause here averaging over the whole volume of colliding 
system is carried out rather than over some fixed central 
region. Note that local peak values may exceed these 
average values by a factor of about 2. 

To get notion of the extent of thermalization achieved 
in the collision process, in Fig. El we show the time 
evolution of the entropy per baryon number (the latter 
is conserved quantity) in central collisions. These cal¬ 
culations were performed without freeze-out in order to 
keep account of the total entropy. We have taken into 
account the entropy of only formed fireball fluid, since 
for the unformed part even its definition is ambiguous. 
Therefore, we deal with an open system, when not all 
the fireball matter is formed. At the early stage of the 
collision (till the time instant of complete overlap of nu¬ 
clei) the entropy quickly rises, which is the evidence of 
nonequilibrium in the system. After that the total en¬ 
tropy either flatten (at E^h < 40A GeV) or even slightly 
drops down (at E\sh > 40A GeV). This entropy flattening 
naturally means thermalization. The temporal dropping 
down of the entropy is related to still continuing produc¬ 
tion of unformed fireball matter. It starts to form only 
later, when the entropy starts to rise after the tempo¬ 
ral fall. Therefore, at the early stage of the reaction we 
deal with an open system. The entropy production in 
the baryon-rich fluids is proportional to the term 

R = (wp • Ut) {Dp — De) — {Dp -\- Dp ), 

if the fireball fluid is still unformed, cf. Eqs. 0, o 
and y^. This result follows from the standard deriva¬ 
tion |96ll of hydrodynamic equation for the entropy flow. 
Therefore, the fall or rise of the entropy of the baryon-rich 
fluids depend on the sign of this R. If i? > 0, the entropy 
rises. This occurs either at Dp = 0, since (up • Ut) > 1, 
or when {up • Ut) S> 1, as it is the case at the initial in¬ 
terpenetration stage. In the region, where the entropy 
drops down, Dp is still nonzero and {up ■ Ut) is not high 
enough, therefore R < 0. In physical terms, the entropy 
drops down, when the energy-momentum radiation into 
the fireball fluid oc — [(wp • wt) -f 1] Dp dominates over 
the internal heating oc [(wp • Ut) — 1] Dp of the baryon- 
rich subsystem. 

As seen from Fig. El at -E'lab < 40A GeV the entropy 
is approximately conserved already after the colliding nu¬ 
clei overlap. These overlaps correspond to peak values in 
Fig. II71 and turning points of the trajectories in Fig. El 
This fact suggests that thermalization in these collisions 
occurs comparatively early. At higher incident energies, 
E\ah > 40A GeV, the complete thermalization happens 
at the comparatively late stage (t > 4 fm/c), when the 
fireball fluid gets formed and the entropy approaches the 
plateau. However, the entropy of the baryon-rich fluids, 
displayed by dashed lines in Fig. El slowly changes al¬ 
ready after the complete overlap of nuclei. It remains 
constant within 10%. This baryon-rich subsystem can 
be considered as an approximately thermalized fluid still 
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interacting with a “bath” of the fireball matter. There¬ 
fore, this overlap time can be approximately taken as a 
time of equilibration in the baryon-rich subsystem at all 
considered here incident energies. 



<fi^> [MeV] 


FIG. 19: (Color online) Dynamical trajectories in the (T, /rs) 
plane for central Pb-|-Pb collisions {h = 2.5 fm) at various 
incident energies. Numbers near the dynamical trajectories 
indicate the evolution time instants in the c.m. frame of 
the colliding nuclei. Bold parts of trajectories are related 
to approximately thermalized baryon-rich subsystem, while 
the thin ones - to yet nonequilibrium evolution. The light- 
grey shaded region corresponds to the boundary of the phase 
transition from the hadronic phase to the QGP, as it was 
estimated in Ref. 1^ . Dotted line is the “experimental” 
freeze-out curve fitted to observed multiplicities in the ap¬ 
proximation of the ideal gas model under condition that 
the energy per hadron is 1 GeV. The star-symbol is the critical 
end-point calculated in Ref. |^. 


Based on the EoS used in simulations, the ns and 
e densities can be recalculated in the temperature- 
chemical-potential representation and displayed 

in the conventional way as dynamic trajectories in the 
plane, see Fig. ^1 In view of the above dis¬ 
cussion, we display approximately equilibrium (for the 
baryon-rich subsystem) parts of the trajectories by bold 
lines, while those yet nonequilibrium - by thin lines, just 
in order to keep track of the system evolution. As seen, 
the compression stage is rather short. For example, at 
Flab = 158A GeV the maximal values of (T) and (itb) 
are reached during the time of about 1 fm/c. After that 
a comparatively fast expansion, accompanied by rather 
slow cooling, starts. At the incident energy of « 30A 
GeV, the ^namical trajectory passes near the critical 
end-point [g , where the order of the deconfinement phase 
transition is changed. In fact, the trajectories in the 
(T,^b) plane strongly depend on the used EoS, even if 
they originate from the same (ub, e) initial state. There¬ 
fore, this should be taken only as a rough estimate of 
the relevant incident energy, since the EoS of the hydro 




FIG. 20: (Color online) Invariant 4-volume corresponding to 
conditions on (solid lines) and on (long-dashed 

lines) for central Pb-|-Pb collisions as function of the incident 
energy. The dotted line on the lower panel is the 4-volume 
At ■ ttR^ ■ 2Rljcm with R — 4: fm, At = 3 fm/c (see the text). 


simulations has not involved any QGP phase transition 
at all. The boundary of the phase transition from the 
hadronic phase to the QGP, estimated in Ref. (light- 
grey shaded region Fig. CH), is displayed just to remind 
on what may happen in this thermodynamic domain, if 
a EoS involves the phase transition. However, this tra¬ 
jectory representation is not quite informative due to the 
following reason. In addition to the uncertainties of map¬ 
ping of nonequilibrium on equilibrium above discussed, 
the displayed time instants correspond to the cm frame 
of colliding nuclei and hence are noninvariant. Therefore, 
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Fig. El does not tell us (in invariant way) how long the 
dense matter survives and in which volume. 

The way to overcome the above difficulties was pro¬ 
posed in [ij. It consists in calculation of an invariant 
4-volume V4 in which a quantity q exceeds a given value 

Q 

V4{Q)= j d^xQ{q-Q). (48) 

This quantity provides a Lorentz invariant measure of the 
space-time region, where the quantity q keeps high value 
q > Q. In Fig. E) this 4-volume is shown for two cases: 
q = and q = \ The first case corresponds 

to the 4-volume summed only over those regions, where 
full stopping (i.e. i? ~ 0 in terms of Eq. (ESJ) has oc¬ 
curred and thermalized ns exceeds certain value, while in 
the second case the 4-volume is summed over all regions 
(i.e. similarly to that in Eq. Ml. including those where 
stopping has not occurred. As seen, for production of the 
matter with comparatively high densities (jib > 4 no ), 
there are certain preferable incident energies, which al¬ 
low to attain the largest 4-volume. In particular, incident 
energies (10-40)A GeV, planned at the new GSI facil¬ 
ity, are favorable for production of equilibrated matter 
with baryon densities higher than 6no. This conjecture 
is heavily based on the hadronic EoS used. However, 
we expect that it is not too far from the truth, because 
we fairly well reproduced the observed stopping power, 
which is of prime relevance to the achieved compression. 
To get an impression of the scale of the corresponding 
4-volume, we draw a short-dashed line presenting the 4- 
volume At-nR'^ ■2R/"fcm composed of Lorentz-contracted 
cylindrical space volume of radius R = Aim. and time in¬ 
terval At = 3 fm/c, here 7cm is the 7-factor of colliding 
nuclei in their c.m. frame. It is of interest that the fall-off 
of the invariant 14 is defined by this 7-factor. 

V. DISCUSSION AND CONCLUSIONS 

In this paper we have presented an extension of the 
3-fluid model for simulating heavy-ion collisions in the 
range of incident energies between few and about 200A 
GeV. In addition to two baryon-rich fluids, which con¬ 
stitute the 2-fluid model, the baryon-free fireball fluid of 
newly produced particles with delayed evolution is incor¬ 
porated. This delay is governed by a formation time, 
during which the fireball fluid neither thermalizes nor 
interacts with the baryon-rich fluids. After the forma¬ 
tion, it thermalizes and comes into interaction with the 
baryon-rich fluids. This interaction is estimated from el¬ 
ementary pion-nucleon cross-sections. 

The hydrodynamic treatment of heavy-ion collisions 
is an alternative to kinetic simulations. The hydrody¬ 
namic approach has certain advantages and disadvan¬ 
tages. Lacking the microscopic feature of kinetic sim¬ 
ulations, it overcomes their basic assumption, i.e. the as¬ 
sumption of binary collisions, which is quite unrealistic in 


dense matter. It directly addresses the nuclear EoS that 
is of prime interest in heavy-ion research. Naturally, we 
have to pay for these pleasant features of hydrodynamics: 
the treatment assumes that the non-equilibrium stage of 
the collision can be described by the 3-fluid approxima¬ 
tion. 

Taking advantage of the modern computers, substan¬ 
tial work has been also done on improvement of numerics 
of the model. In particular, it was found that the numer¬ 
ical diffusion of the computation scheme should be as 
much as possible avoided in order to get reliable results 
in hydro simulations. To avoid the numerical diffusion 
the parameters of the scheme should be carefully opti¬ 
mized, as it is described in App. El 

The main unknowns of the present simulations can be 
briefly summarized as follows: the equation of state and 
“cross sections”. The EoS is an external input to the 
calculation and thus can be varied. Our goal is to find 
an EoS which in the best way reproduces the largest body 
of available observables. The “cross sections” are equally 
important. They determine friction forces between fluids 
and hence the nuclear stopping power. At present we 
have at our disposal only a rough estimate of the friction 
forces (cf. Ref. for the friction of baryon-rich fluids 
and Sect. mn for the friction with the fireball). There¬ 
fore, we have to fit the friction forces to the stopping 
power observed in proton rapidity distributions. 

In this paper we use a purely hadronic scenario in our 
simulations. This means a very simple EoS con¬ 
structed under the requirement that it reproduces satu¬ 
ration of the cold nuclear matter. This EoS is a natu¬ 
ral reference point for any other more elaborate EoS. We 
have found out that the original friction between baryon- 
rich fluids, incorporated in the model accordingly to es¬ 
timate of Ref. , is evidently insufficient to reproduce 
observable stopping of nuclear matter. Therefore, the 
original friction was enhanced by means of the tuning 
factor fh, cf. Eq. lEJ, which is plotted in Fig. |21 As 
seen, the enhancement factor turned out to be unexpect¬ 
edly large. However, it complies with earlier results of the 
2-fluid model with mean mesonic fields . There it was 
found out that for the proper reproduction of the data 
on heavy-ion collisions in the energy range from 0.4A 
to 0.8A GeV, where the purely hadronic scenario is cer¬ 
tainly applicable, the enhancement by a factor of 3 was 
required. From this large enhancement we infer that the 
rough estimate of the friction [24] is really too rough and 
that further more elaborate microscopic calculations of 
the friction are required. In doing this at lower energies, 
we have to rely on pure hadronic effects, e.g. on the con¬ 
cept of the hadron liquid rather than hadron gas. At 
higher incident energies the friction enhancement could 
be associated with an indirect manifestation of the onset 
of the phase transition, alt hough the used EoS does not 
involve it directly. In Refs. |lO(lllo4ll02| it was argued 
that just above the critical temperature, the quark-gluon 
system is still a strongly interacting matter rather than 
a gas of perturbative partons. This strong interaction 
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anomalously enhances interaction cross sections in this 
phase. 

With this simple hadronic EoS we have succeeded to 
reasonably reproduce a great body of experimental data 
in the incident energy range Eiab — (1-160)^ GeV. The 
list includes proton and pion rapidity distributions, pro¬ 
ton transverse-mass spectra, rapidity distributions of A 
and A hyperons and antiprotons, elliptic flow of protons 
and pions (with the exception of proton V 2 at 40A GeV), 
multiplicities of pions, A and A hyperons, AT"*" mesons 
and (j) mesons, 5^ and multistrange hyperons. How¬ 
ever, we have also found out certain problems. The calcu¬ 
lated yield of K~ is approximately a factor of 1.5 higher 
than that in the experiment. We have also failed to de¬ 
scribe directed transverse flow of protons and pion at 
Eiah > 40A GeV. This especially concerns pions, vi of 
which experimentally anticorrelate with that of protons, 
while it correlates in our calculations. This failure appar¬ 
ently indicates that the used EoS is too hard and thereby 
leaves room for a phase transition, which always makes 
the EoS softer. 

The used hadronic EoS allows modifications relating 
to different incompressibilities (K) of the nuclear mat¬ 
ter. We performed simulations with soft {K = 210 and 
235 MeV) and hard {K = 380 MeV) EoS’s. We have 
found out that the data certainly prefer soft EoS. All the 
calculations presented in this paper correspond to K = 
210 MeV. 

The analysis of the global evolution of central Pb-f Pb 
collisions with present hadronic scenario shows in partic¬ 
ular that incident energies (10~40)A GeV, planned at the 
new GSI facility, are favorable for production of equili¬ 
brated matter with baryon densities 6 times higher than 
the normal nuclear density. This conclusion is heavily 
based on the simple hadronic EoS used in our calcula¬ 
tions. However, we expect that it is not too far from 
the truth, because we fairly well reproduced the ob¬ 
served stopping power, which is of prime relevance to the 
achieved compression. Moreover, in this energy range 
dynamical (T,/is) trajectories of nuclear collisions pass 
in the vicinity of the critical end-point, as it is estimated 
in Ref. |^. Note that the presently used EoS does not 
involve any phase transition. In view of the latter, it is 
intriguing that discrepancies of the present hadronic sce¬ 
nario with the data (i.e. in the form of the A rapidity 
distribution, proton and pion vi and V 2 , and strangeness 
production) are most clearly seen at incident energies 
about 40A GeV. 

Our present experience shows that in order to conclude 
on the relevance of a particular EoS the whole available 
set of data should be analyzed in a wide incident energy 
range with the same fixed parameters of the model. This 
is one of the main conclusions of our paper. Indeed, any 
particular piece of data can be apparently fitted by means 
of fine tuning of these parameters, while to fit the whole 
set of them, one needs a really good EoS. In the present 
paper we analyzed various data only fragmentarily, just 
to fit the parameters and to give an overview of the re¬ 


sulting predictions, since the framework of a single paper 
does not allow us to do more. Our preliminary conclu¬ 
sion is that the considered hadronic EoS is not perfect at 
higher incident energies Ei^h > 20A GeV because of the 
above mentioned problems. Detailed analysis of particu¬ 
lar pieces of available data will be reported in forthcom¬ 
ing papers. We are also going to extend our analysis to 
other EoS’s involving phase transition to the QGP. 

Web page 

The source code of the 3-fluid model is publically avail¬ 
able on the Web page http://theory.gsi.de/~mfd/. Free 
downloads of the code are available from there. Brief in¬ 
structions on how to run this code are also given on this 
page. 
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APPENDIX A: NUMERICS 

The relativistic 3D code for the above described 3-fluid 
model was constructed by means of modifying the ex¬ 
isting 2-fluid 3D code of Refs. mill. The numeric 
scheme of the code is based on the modified particle-in¬ 
cell method |4^ l4^ , which is an extension of the scheme 
first applied in Los-Alamos ^3- Taking advantage of 
modern computers, we have comprehensively tested this 
numeric scheme in order to find optimal parameters of 
the calculation. 

In the particle-in-cell method the matter is represented 
by an ensemble of Lagrangian test particles. They are 
used for calculation of the drift transfer of the baryonic 
charge, energy and momentum. In the present scheme 
the test particle has a size of the cell. Therefore, when a 
single test particle is moved on the grid, it changes quan¬ 
tities in 8 cells, with which it overlaps. These spatially 
extended particles make the scheme smoother and hence 
more stable. The transfer due to pressure gradients, fric¬ 
tion between fluids and production of the fireball fluid is 
computed on the fixed grid (so called Euler step of the 
scheme). Simulation is performed in the frame of equal 
velocities of colliding nuclei. In view of this combined 
nature of the numerical scheme there are parameters of 
the test particles and the grid, which are summarized in 
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the Table im 

The basic scale of the numerical scheme is chosen along 
the X direction (the transverse direction in the reaction 
plane). According to the preset number of cells in the x 
direction, L^, the maximal impact parameter &niax, which 
is foreseen for the simulation, and number cells free of the 
matter, reserved for transverse expansion, the Ax 

step of the grid is determined as follows 

Aa; = max Up + i?t + fomax} /— L^“^®®), 

where Up and Rt are r adii of the projectile and target 
nuclei, respectivelyAnother transverse step Ay (in 
the out-of-reaction-plane direction) is taken equal to Ax: 
Ay = Ax. Az and At are determined by means of preset 
ratios Ax/At and Ax/Az. The preset numbers of cells in 
the y and z directions {Ly and L^) are used to finally de¬ 
fine the size of the grid. The simulation takes place only 
in the semi-space y > 0 because of the symmetry of the 
system with respect to the reaction plane. The number 
of baryon-rich test particles {N//// = A^p“jectiie + {Vtarget) 
is specified at the initialization step of the simulation. As 
for produced baryon-free fluid, the actual total number 
of corresponding test particles depends on the incident 
energy, the size of the cell and the number of produced 
fireball test particles per cell and per time step (Nf) and 
is determined only upon the completion of the simulation. 
The maximal total number of baryon-free test particles 
iV^ot® should be larger than the actual total number. 
is used to reserve memory for baryon-free test particles. 


Slab, A GeV 

1 T 10 

20 -h 158* 

158** 

Lx 

150 

320 

320 

Ly 

60 

130 

130 

L, 

240 

820 

820 

T free 

5 

43 

43 

Ax/At 

3.5 

3.5 

4 

Ax/Az 

1 

1 

1.2 


2 • 10® 

15 ■ 10® 

15 • 10® 

Nf 

2 

2 

2 


3-10® 

25 ■ 10® 

25 • 10® 


TABLE II: Basic parameters of the numeric scheme used for 
simulations of Au-|-Au and Pb-|-Pb collisions at various inci¬ 
dent energies Aiab. 


In order to get reasonable accuracy in the simulation, 
the following basic requirements should be met: 

(i) The grid should be extended enough to prevent escap¬ 
ing the matter beyond this region before it gets frozen 
out. 

(ii) The grid in the beam {z) direction should be fine 
enough for a reasonable description of the matter of ini¬ 
tially Lorentz-contracted nuclei. From the practical point 
of view, it is desirable to have more than 30 cells on the 
Lorentz-contracted nuclear diameter. 

(iii) The well-known Courant-Friedrichs-Lewy criterion 


states that the ratios of the space-grid steps to the time 
step (e.g.. Ax/At) should be larger than 1 in order to 
have a consistent and stable algorithm for solving hy¬ 
perbolic partial-differential equations. To avoid numeri¬ 
cal diffusion, this ratios should be taken optimal. As it 
was found in 1 -dimensional simulations of exactly solv¬ 
able problems the optimal range of these ratios is 
2.5 < Ax/At < 6 with the preferable Ax/At ~ 3.5, 
minimizing the numerical diffusion. This fact dictates 
the choice of equal-step grid in all directions (Ax : Ay : 
Az = 1 : 1 : 1), in spite of Lorentz-contraction of inci¬ 
dent nuclei, which is quite strong at high energies. This 
choice makes the scheme isotropic with respect to the 
numerical diffusion. However, in view of (i) and (ii) re¬ 
quirements it makes the grid too fine in the transverse 
directions and thus results in high memory consumption. 
The need of the equal-step grid in all directions for rela¬ 
tivistic hydrodynamic computations within conventional 
1-fluid model was pointed out in Ref. ^3. As it was 
demontrated there, the matter transport becomes even 
acausal if this condition is strongly violated. 

(iv) The number of test particles per cell should be 
high enough to avoid large numerical fluctuations in drift 
transfer. In practice, this means that this number should 
be not less than 3 during the whole evolution of expand¬ 
ing matter, till its complete freeze-out. 

The final check of these requirements is possible 
only upon completion of the simulation. The Table QH 
presents the optimized set of parameters for calculations 
of Au-|-Au and Pb-|-Pb collisions at various incident en¬ 
ergies. This set was determined by multiple test runs 
of the code. The reaction Pb -I- Pb at i?iab = 158A 
GeV is the most memory and time consuming case con¬ 
sidered in simulations. To completely fulfill the require¬ 
ment (iii) in this calculation we need 8 GB of memory 
for the central collision (see column 20 -L 158* in Table 
ini). In order to comply with memory restrictions, we still 
had to take slightly unequal steps in different directions: 
Ax : Ay : Az = 1.2 : 1.2 : 1, Ax/At = Ay/At = 4 
and Az/At = 3.3 (see column 158** in Table in)l . for 
noncentral Pb -I- Pb collisions at Aiab = 158A GeV. 

The freeze-out criterion is checked at the La- 
grangian step of the scheme. It is checked in the cells, 
which overlap with the considered test particle. As men¬ 
tioned above, the test particle has a volume of the cell. 
Therefore, each test particle overlap with 8 cells. If the 
freeze-out criterion is met in all these cells, this test parti¬ 
cle is frozen out and removed from further hydrodynamic 
evolution. Thus, the frozen-out test particles are those 
droplets mentioned in Eq. for the spectrum of ob¬ 
servable hadrons. Precisely due to this extended chec king 
in 8 cells, the actual freeze-out energy density [TTH 
participating in the freeze-out criterion (EDI differs from 
that used in the code £f°/®. According to our estimate, 
the actual freeze-out energy density e^z is approximately 
twice as lower than ef^J^®. This estimate was checked nu¬ 
merically. The reason behind this extended checking in 
8 cells is as follows. Had we checked the criterion only 
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in a single cell, where the test particle is located, we 
would prevent any hydrodynamic expansion of the sys¬ 
tem. Indeed, if the matter, during its expansion, starts 
to fill in a cell, which was empty before that, it occupies 
less than 1/3 of this cell during the single time step, be¬ 
cause the time step should be less than the space step 
(cf. Ax/At = 3.5) for this numerical scheme. Then, all 
the densities in this newly filled-in cell are quite low at 
this time step. Therefore, the freeze-out criterion is cer¬ 
tainly fulfilled in this cell. Had we confined ourselves to 
checking only this cell, we would immediately freeze the 
matter in this cell out, thus preventing hydrodynamic 

I 


expansion of the matter. 

In order to avoid formation of bubbles of frozen out 
matter inside the dense environment, we introduced ad¬ 
ditional criterion of the freeze-out: at least one of the 
above discus sed c ells should be empty, i.e. contain no 
test particles This additional criterion means that 

the analyzed test particle is located near the surface of 
the system, provided, of course, that the number of test 
particles per cell is large enough (cf. condition (iv)), 
which excludes origination of such empty cells due to 
fluctuations. 


APPENDIX B: RESONANCE DECAYS 


1. General Consideration 


Contribut ion o f resonance decays into stable particle spectra has been studied long ago |103| and later analyzed in 
more detail |104| . To m ake the p aper self-contained, we briefly present here formulas for these decays, which in fact 
are the same as those in but differ from those in presentation, which, from our point of view, is advantageous 

for numerical realizations. 

Consider a hadronic resonance R with degeneracy factor gn and a spectral function A^(s), where s is the resonance 
mass squared, which decays into n particles 

R —> 1 -p 2 -j- ... -p n 


of masses mi, m 2 ,..., m„, through a Jth channel of its decay with a branching ratio bj. Then the distribution of the 
”1” particle, resulting from this decay, is as follows 


El- 1 


d^Pi 


= bj 




ds Aii{s) 


d^PR (^^d?Nti{s,pR) 


Er 


d^PR 


El 


d?v[^~ 


■1+2+...+n) 


{s,Pr,Pi) 


d^Pi 


Here, the distribution of the R resonance, produced by the hydrodynamic computation, is 

d^NR{s,pR) prV u-pR 


Er- 


d^PR 


(27r)3 exp[/3(u • pr - pr)] ± 1 ’ 


(Bl) 


(B2) 


where P is a small proper (i.e. in the rest frame) volume of the fluid element, gj^ is degeneracy of the R resonance, 
u is the hydrodynamic 4-velocity, j3 is the inverse temperature, pR is the chemical potential, pr is the 4-momentum 
of the R resonance, s = the upper (lower) sign in this expression corresponds to baryonic (mesonic) resonances. 
The distribution of the ” 1” particle, resulting from decay of a single R resonance, is expressed as follows 


El 


+r(^^i+^+-+")(s,Pi,,Pi) 

d^Pi 



(B3) 


unde r assumption that the matrix element of the decay is constant all over the available n-particle phase-space volume 
|l05| . which, in its turn, reads 


Rn{s) 



(B4) 


In order to obtain the total contribution of the R decay into ” 1” particle spectrum, we should sum over all J branches 
of the R decay, in which the ”1” particle appears, and take into account that particles identical to ”1” can also be 
among residue particles 2, 3,..., n 
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d?pi 
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(B5) 
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Here Ni is the multiplicity of ” 1” 


particles in the Jth channel of the R decay. Here and below, Eji = (s + p|j) 


1/2 


and Ei = (mf + 

Fortunately, the the major part of integrations involved in Eqs. m and {HU) can be performed analytically. We 
explicitly present results only for two-particle and three-particle decays, which are dominant for resonances under 
consideration. 


2. Two-Particle Decays 


The result for a two-particle decay of a i? resonance is as follows 


E, 


(Epi 


= b. 


gnV 




1 _ 

— J (mi+m 2 )' 


A / \ OO 


n—1 


lEl 


n(3ER + l f ” 

f 0^2 Wfl - g-R)} 

Mr 


where the upper (lower) sign in this expression corresponds to baryonic (mesonic) resonances. 


^ 2 (s) = 

2s 


is the 2-particle phase-space volume and 


E^ = 
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« 2 m 2 


(s -I- mM RT-l) {r ■ Pi) 4= ml,ml) [(w • pi)^ — mf] 


21 1 / 2 ' 


Here and below 

Hx, y,z) = [x - y - z)'^ - Ayz 

is the standard kinematic function |10 !t| . To arrive to this result, the following expansion 
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exp[P{u-pR - Hr)] ± 1 


^ (=Fl)" ^ exp {-nP{u ■ pr - pr)} 
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(BIO) 


was used, the advantage of which is that it is rapidly convergent, when u-pR > pr. The opposite limit, i.e. u-pR < pr, 
corresponds to low temperatures and hence hardly contributes to the resonance production. The n = 1 term in this 
expansion, as well in Eq. (IB6II . corresponds to the classical Jiittner distribution. 


3. Three-Particle Decays 


The result for a three-particle decay of a i? resonance is as follows 
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where the upper (lower) sign in this expression again corresponds to baryonic (mesonic) resonances, 

= min ^ , X = max ^X ; 

Xmax = , X^zn = {^2 + TTlsf , 

= s + mi-2ERiu-pi)±2{E]i-[{u-pif 

A* = A {X^,ml,ml) , 

Rsis) = (I) iy {s,X,ml) {X,ml,rnl) , 

9{X'^), 9{X~), etc. are step functions which define the accessible kinematic region. 


(B12) 

(B13) 

(B14) 

(B15) 
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